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ABSTRACT 

A  boundary  integral  formulation  for  the  analysis  of 
circular  plate  bending  under  lateral  loads  is  developed  using 
Green’s  functions.  The  formulation  specifically  applies  to 
annular  plates  with  arbitrary  boundary  conditions.  The  plate 
bending  solution  in  the  plastic  range  is  determined  using  a 
numerical  method  of  incremental  loading.  A  computer  program  to 
perform  the  required  calculations  was  developed  and  is 
presented.  Results  for  three  case  studies  are  included  and 
compared  with  results  obtained  by  other  methods.  Plate 
behavior  in  the  elastic  range  is  in  excellent  agreement  with 
other  analytical  sol'utions,  and  in  the  plastic  range  is  in 
reasonable  agreement  with  published  results  obtained  using  a 
finite  element  method. 


Thesis  Supervisor:  Dr.  Amiram  Moshaiov 

Assistant  Professor  of  Ocean  Engineering 


/\f#l/fr/  P0ST&#&2>u#r*  SC-Atsf 
/tfSJL AT*  fey  (2#  43  9Y2  -  d 

a/pc>  <22$  - 

r 


ACKNOWLEDGEMENTS 


I  wish  to  express  my  deepest  appreciation  for  the  help  and 
support  that  I  received  from  my  thesis  advisor,  Amiram 
Moshaiov,  throughout  every  stage  of  development  of  thi3  thesis. 
His  ability  to  explain  complex  concepts  very  simply,  his 
enthusiasm,  his  unending  patience  and  his  willingness  to  set 
aside  work  at  a  moments  notice  to  lend  assistance  to  me  as  well 
as  other  students,  all  earn  him  my  highest  marks  as  a  teacher. 
The  countless  hours  he  spent  with  me  have  contributed  more  to 
my  education  than  any  other  experience  at  MIT.  For  the 
pleasure  of  working  under  his  supervision,  I  will  be  forever 
grateful . 

To  Professor  Jerome  Connor  of  the  Civil  Engineering 
Department,  I  must  offer  thanks  for  sparking  my  interest  in  the 
this  general  thesis  area  as  I  was  struggling  through  his  course 
in  Finite  and  Boundary  Element  Methods.  His  enthusiasm  in  and 
knowledge  about  these  topics  have  been  an  inspiration  to  me. 

Finally,  I  wish  to  express  my  gratitude  to  CAPT  Clark 
Graham,  USN  and  CDR  Dave  Whiddon,  USN  who  administered  the 


3 


X I  1 1 A  program  during  my  final  two  years  at  MIT.  Their 
understanding,  patience,  support  and  positive  approach  made  my 
educational  experience  here  a  profitable  and  enjoyable  one. 


4 


TABLE  OF  CONTENTS 


Page 


1  INTRODUCTION .  12 

2  DEVELOPMENT  OF  THE  GOVERNING  DIFFERENTIAL  EQUATION .  15 

2.1  Moment-Curvature  Relationships .  15 

2.2  Equilibrium  Equations .  21 

2.3  Differential  Equation  for  Circular  Plate  Bending....  23 

3  A  BOUNDARY  INTEGRAL  METHOD .  25 

3.  1  General  Formulation .  25 

3.2  Beam  Bending  Problem .  27 

4  DEVELOPMENT  OF  THE  FORMULATION  FOR  CIRCULAR  PLATES .  32 

4.1  Description  of  Plate  Geometry .  3  2 

4.2  Development  of  Integral  Equation .  3  4 

4.3  Selection  of  a  Green's  Function .  3‘J 

4.4  Solution  of  Integral  Equation .  40 


A 


5 


TABLE  OF  CONTENTS  fccntd.  ; 

Page 

5  NON-LINEAR  BEHAVIOR .  45 

5.1  Plastic  Stress-Strain  Relationships .  45 

5.2  Plastic  Moment-Curvature  Relationships .  46 

5.3  Equilibrium  Equation  With  Plasticity .  48 

5.4  Boundary  Integral  Formulation  With  Plasticity .  49 

6  NUMERICAL  SOLUTION .  58 

6.1  Incremental  Load  Method .  58 

6.2  Computer  Program  Description .  62 

6.3  Numerical  Methods .  65 

7  RESULTS  AND  DISCUSSION .  66 

7.1  Simply  Supported  E 1  as t o- P 1  as t i c  Annular  Plate .  66 

7.2  Clamped  E  1  as  t  o  -  P  1  as  t  i  c  Annular  Plate .  69 

7.3  Simply  S uppo r t ed / Gu i ded  E  1  as t o-P  l  as t  i  c 

Annular  Plate .  72 

7.4  Discussion  of'  Results .  77 

8  .SUMMARY  AND  CONCLUSIONS 

8.1  Summarv .  79 

8.2  Conclusions .  79 

8.3  Recommendations  for  Future  Work .  80 


6 


TABLE  OF  CONTENTS  (eontd) 

Page 


Appendix  A.  Selected  Green's  Functions .  82 

Appendix  B.  Computer  Programs .  94 

List  of  References .  133 


7 


LIST  OF  FIGURES 

Figure  Page 

2-1  Deflection,  slope  and  curvature  relations  hips 

in  circular  plate  bending .  IV 

2-2  Differential  circular  plate  element .  19 

2- 3  Shear  forces  and  bending  moments  acting  on 

differential  plate  element .  22 

3- 1  Examples  of  Betti’s  reciprocal  work  theorem .  29 

4- 1  Symmetrically  loaded  annular  plate .  3  3 

4-2  Plate  geometry  for  the  selected  Green's  function....  40 

6-1  Uniaxial  stress -strain  curve  for  elastie- 

perfectly  plastic  material .  62 

6- 2  Computer  program  flow  diagram .  64 

7- 1  Simply  supported  annular  plate  problem 

desc  r  lption .  67 

7-2  Deflection  of  simply  supported  annular  plate 

at  the  onset  of  yielding .  68 

7-3  Deflection  of  simply  supported  annular  plate 

in  the  plastic  range .  68 

7-4  Plastic  zones  in  simply  supported  annular 

plate .  69 

7-5  Clamped  annular  plate  problem  description .  70 

7-6  Deflection  of  clamped  annular  plate  at  the 

onset  of  yielding .  71 

7-7  Deflection  of  clamped  annular  plate  in  the 

plastic  range .  71 


8 


LIST  OF  F IGURES ( con t  d ) 


Figure  Page 

7-8  Plastic  zones  in  clamped  annular  plate .  72 

7-9  Simply  s uppc r t ed/ gu i ded  annular  plate 

problem  description .  73 

7-10  Deflection  of  simply  s uppo r t ed/ gu i d ed  annular 

plate  at  the  onset  of  yielding .  74 

7-11  Deflection  of  simply  s uppor t ed/gu i ded  annular 

plate  in  the  plastic  range .  74 

7-12  Plastic  zones  in  simply  suppo r t ed/ gu i ded 

annular  plate .  75 

7-13  Simply  supported  circular  plate  problem 

description  from  Moshaiov  and  Vorus .  75 

7-14  Deflection  of  simply  supported  circular  plate 

in  the  plastic  range  from  Moshaiov  and  Vorus .  76 

7-15  Plastic  zones  in  simply  supported  circular 

plate  from  Moshaiov  and  Vorus .  76 


9 


b 

c 

D 

d 

E 

h 

I 

K  r 
Kfl 

L 

M  ,  MG 

*r.MGr  - 

Mrp 

Mr 

^a.^Ge 

MflP 

q,qg 

Qr.0Gr  “ 

Or 

0 . 9G 


NOMENCLATURE 


plate  outer  radius 
plate  inner  radius 

plate  outer  radius  for  Green’s  function 
flexural  rigidity 

plate  inner  radius  for  Green’s  function 
Young’s  modulus 
plate  thickness 

area  moment  of  inertia  of  beam  cross  section 
radial  curvature 
tangential  curvature 
beam  length 

elastic  moments  for  beam  bending  problem 
elastic  radial  moments  per  unit  length 
plastic  radial  moment  per  unit  length 
actual  radial  moment  per  unit  length 
elastic  tangential  moments  per  unit  length 
plastic  tangential  moment  per  unit  length 
shears  for  beam  bending  problem 
elastic  radial  shears  per  unit  length 
actual  radial  shear  per  unit  length 
distributed  loads  per  unit  area 
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radial  radius  of  curvature 
tangential  radius  of  curvature 
radial  distance  from  plate  center 
radial  position  of  ring  load 
vield  stress 
deflections 
distance  along  beam 

lateral  distance  from  plate  center 

plastic  strain  increment 

total  radial  strain 

elastic  radial  strain 

plastic  radial  strain 

total  tangential  strain 

elastic  tangential  strain 

plastic  tangential  strain 

yield  strain 

s 1  opes 

plate  boundary 

angular  position  in  tangential  direction 
biharmonic  operator- 

loads  per  unit  length  for  beam  bending  problem 

magnitude  of  Green’s  function  ring  load 

equivalent  stress 

radial  stress 

tangential  stress 

Poisson's  ratio 

plate  domain 
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CHAPTER  1 


INTRODUCTION 


The  solution  of  plate  bending  problems  in  the  plastic 
range  has  been  approached  using  finite  element,  and  more 
recently,  boundary  integral  methods  (also  called  boundary 
element  methods).  Unlike  the  finite  element  method,  which 
discretizes  the  inside  of  the  plate  into  a  number  of  small 
elements,  the  boundary  integral  method  discretizes  only  the 
plate  boundary,  thereby  reducing  the  order  of  the  problem  by 
one . 


Many  different  approaches  have  been  used  to  formulate 
plate  bending  problems  using  boundary  integrals,  with  perhaps 
the  most  attractive  proposed  by  Stern''^  and  8ezine<-2>,  These 
authors,  working  independently,  developed  a  direct  boundary 
integal  method  using  Green's  functions  to  solve  the  general 
plate  bending  problem  in  the  elastic  range.  This  work  has  been 
extended  by  Moshaiov  and  Vorus^^  to  include  plastic  behavior. 

Symmetry  is  commonly  used  to  reduce  the  size  of  the  system 
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of  equations  needed  to  analyze  a  symmetric  structure.  In  the 


case  of  ax i symme t r i ca 1 1 y  loaded  circular  plates  the  problem 
becomes  one-dimensional.  For  the  linear  elastic  case,  it  has  a 
closed  form  analytic  solution.  However,  some  authors  have 
instead  treated  this  circular  plate  problem  as  a  two- 
dimensional  problem,  perhaps  for  the  sake  of  demonstration 
(Kamiya  and  Sawki*-^^).  They  have  used  symmetry  to  allow 
discretization  of  a  portion  of  the  boundary  instead  of  the 
entire  boundary. 

This  thesis  develops  a  general  formulation  for  solving 
annular  plate  bending  problems  using  boundary  integrals  that 
reduces  the  problem  to  essentially  a  one- d imens i ona 1  problem  by 
the  use  of  a  "ring"  type  Green’s  function.  Such  a  formulation 
is  most  useful,  inspite  of  the  closed  form  solution,  as  it 
allows  the  analysis  of  plates  with  different  boundary 
conditions  and  loads  with  one  algorithm.  Moreover,  it  permits 
a  solution  in  the  plastic  range,  which  is  not  possible  with  a 
conventional  analytic  approach.  In  addition,  the  simplicity  of 
the  one-dimensional  formulation  given  here  has  a  unique 
significance  for  the  teaching  of  boundary  element  methods. 

Butterfield’'^  has  developed  a  boundary  element 
formulation  for  the  one-dimensional  elastic  beam  bending 
problem.  Here,  a  similar  approach  is  taken  for  the  annular 
plate.  It  is  also  extended  to  include  non-linear  behavior 
using  an  incremental  load  method  as  outlined  by  Moshaiov  and 


Vorus.  Numerical  results  for  three  different  annular  plate 
configurations  are  presented  and  compared  to  results  obtained 
using  other  methods.  A  discussion  of  these  results  follows, 
well  as  suggestions  for  future  work  in  this  area. 
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CHAPTER  2 

DEVELOPMENT  OF  THE  GOVERNING  DIFFERENTIAL  EQUATION 


This  chapter  reviews  the  derivation  of  the  governing 
differential  equation  for  a  circular  plate  subjected  to  a 
symmetrically  distributed  lateral  load.  The  derivation  is 
based  on  that  provided  by  Timoshenko''*^  ,  using  Kirchoff’s 
assumptions.  It  is,  therefore,  applicable  only  to  the  linear- 
elastic  case.  Treatment  of  non-linear  behavior  is  addressed  in 
Chapter  5 . 


2 . 1  Moment-Curvature  Relationships 

The  first  step  in  developing  the  governing  differential 
equation  for  a  circular  plate  is  to  find  a  relationship  between 
moment  and  curvature.  Curvature  in  the  radial  direction  (Kr) 
and  the  tangential  direction  (Kg)  for  a  symmetrically  loaded 
circular  plate  is  found  using  geometrical  considerations.  Kr, 
which  is  the  inverse  of  the  radius  of  curvature  in  the  radial 
direction  (Hr),  can  vary  as  a  function  of  the  distance  (r)  from 
the  plate  center.  Examining  the  plate  of  Figure  2-1  at  a  small 
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radial  distance  dr  away  from  r,  the  slope  (♦)  can  be  expressed 
in  terms  of  the  radial  distance  r  and  the  deflection  (wi  as 
foil ows 

dw 

♦  = -  (2-1-1) 

dr 

The  radial  radius  of  curvature  is 


giving  an  expression  for  the  radial  curvature  in  terms  of  w 


1 

d  * 

,  2 

d  w 

R 

r 

dr 

dr  2 
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Figure  2-1.  Deflection,  slope,  and  curvature  relatio 
in  circular  plate  bending 


Due  to  symmetry,  the  tangential  radius  of  curvature  is 
constant  at  any  given  radius  r,  and  for  small  deflections  can 
be  expressed  as 


giving  an  expression  for  the  tangential  curvature 
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1 

* 

1  dw 

II 

CD 

r9 

r 

r  dr 

To 

relate  moment 

t  0 

curvature , 

cons i der 

the  small  plate 

e  1  emen  t 

illustrated  in 

Figure  2-2. 

Define  Mr 

and  Mg,  the 

radial 

and  tangential 

moments  per  unit  length 

respectively,  as 

fol lows 

h/2 

M 

r 

o  z  dz 
r 

-h/2 

(2-1-3 

rh/2 

11 

CD 

a  z  dz 

where  h  is  the  plate  thickness.  Assuming  the  plate  is  thin  and 
hence  experiences  a  two-dimensional  stress  state,  Hooke’s  law 
yields  the  following  expressions  relating  stresses  to  strains 


a 

r 


E 

( 1-v 2  ) 


a 


8 


E 

(  1-u  2  ) 


(2-1-4  ) 


where  E  is  Young’s  modulus,  v  is  Poisson’s  ratio,  and  cre  and 
r ge  are  the  elastic  radial  and  tangential  strains,  respec¬ 
tively.  In  Chapter  5,  which  discusses  non-linear  behavior,  a 
distinction  will  be  made  between  the  elastic  strain  and  the 
total  strain.  In  the  plastic  range,  the  total  strain  includes 
elastic  and  plastic  components,  in  accordance  with  the 
following  expressions: 
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.h/2 

"  i  ‘r 

) 

■j  Vs  /  o 


V  £fl  z  dz 


(2-1-5: 


(l-»  ) 


h/2 

18  r  J 

■  _  u  /  n 


Finally,  from  geometrical  considerations,  strain  and  curvature 
are  related  by 


(2-1-6) 


which,  when  substituted  into  Eq.  (2-1-5),  give  the  desired 
moment-curvature  relationship  for  a  circular  plate 


d  w  v  dw 


dr  r  d r 


1  dw 

D - 

r  dr 


where  D,  the  flexural  rigidity,  is  given  by 


(2-1-7) 


12  ( 1 -v  ) 


2-1-8) 
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2 . 2  Equilibrium  Equations 

For  the  momen t -cur va t ur e  relationships  to  be  useful,  an 
equilibrium  equations  in  terms  of  moments  and  shears  must  be 
found.  The  governing  differential  equation  is  then  developed 
by  substituting  the  expressions  relating  moment  to  curvature 
into  the  equilibrium  equation. 


To  find  an  equilibrium  equation,  consider  first  the 
laterally  loaded  circular  plate  element  illustrated  in  Figure 
2-3,  where  Qr  is  the  radial  shearing  force  per  unit  length  and 
q  is  the  distributed  load  (force  per  unit  area).  The  couples 
acting  along  edges  ab  and  cd  due  to  radial  bending  moments  are 


ab  : 

cd : 


d  8 


dr 


r  *  dr 


d  8 


The  shear  forces  acting  along  the  edges  are 


ab :  Q  r  d  8 


cd : 


0  + 
r 


dr 


r  1  dr 


d  8 
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figure  ii-3.  Shear  forces  ana  moments  act  ing  on  Jii  iereuU 

plate  element. 

Because  of  symmetry,  there  are  no  shear  forces  acting 
along  edges  ad  and  be,  leaving  only  the  tangential  bending 
moments 

M0  dr 

These  couples  can  be  resolved  into  components  acting 
perpendicular  to  the  axis  ro  (which  are  equal  and  opposite  and 
thus  cancel),  and  components  acting  along  the  ro  axis  which 
have  a  total  magnitude  of 


Mq  dr  d8 

Summing  up  all  the  couples  and  neglecting  the  small  change  in 
shear  force  across  the  element  (i.e.  dropping  the  d Q  terms) 
gives  the  following  equilibrium  equation 

dM  |  . 

Mr  +  - —  dr  |  |  r  +  dr  |  d8  -  r  d8 

Ma  dr  d9  +  Q  r  d9  dr  =  0 

r 


This  equation  is  further  simplified  by  ignoring  higher  order 
terms  to  give 


M 

r 


M  +  Q  r  =  0 

8  r 


Also,  from  equilibrium 


in  the  z  direction 


(2-2-1  ) 


dQ 


dr 


(  2  -  2  -  2  ) 


2 . 3  Differential  Equation  for  Circular  Plate  Bending 

Having  found  the  desired  equilibrium  equations,  the  moment 
curvature  relationships  of  Section  2.1  are  now  used  to  develop 
the  governing  differential  equation  for  the  circular  plate 
bending  problem.  Specifically,  substituting  Eq.  (2-1-7!  into 
Eq.  12-2-1)  results  in  a  third  order  differential  equation  to 
describe  the  response  of  a  circular  plate  to  lateral  loading 
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d4  2d3  Id2  Id 

7  -  7  4  ~  ?  ~~  2  7  3  d  5  ) 

dr  rdr  rdr  r  d  r 
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CHAPTER  3 

A  BOUNDARY  INTEGRAL  METHOD 


This  chapter  describes  a  boundary  integral  method,  and 
how  it  can  be  used  to  provide,  a  general  method  of  solution  of 
circular  plate  bending  problems.  A  one-dimensional  beam 
bending  problem  is  used  to  illustrate  this  method. 

3 . 1  General  Formulation 

Solution  of  Eq.  { 2 - 3 - 3  )  analytically  as  a  boundary  value 
problem  by  applying  known  boundary  conditions  is  possible 
(Timoshenko).  However,  this  approach  requires  special 
attention  for  each  plate  configuration  and  load  distribution. 
This  difficulty  is  compounded  when  considering  problems 
involving  plastic  behavior  which  should  be  solved  numerically 
due  to  the  non-linear  behavior  in  the  plastic  range. 

Therefore,  a  general  method  is  sought  wherein  the  same 
procedure  can  be  used  to  solve  a  variety  of  problems,  thus 
lending  the  problems  to  computer  solution  techniques.  One  such 
method,  adopted  for  this  work  to  examine  circular  plate 


bending,  is  the  boundary  integral  method. 


Stern  and  Be2ine,  working  independently,  obtained  a 
general  formulation  for  plate  bending  problems  in  terms  of 
boundary  integrals.  These  integrals  involve  displacements, 
slopes,  bending  moments,  and  shears  on  the  plate  boundary. 

Both  authors  used  the  following  reciprocal  identity,  which  can 
be  derived  using  Green’s  identity: 


'  1 

4 

w_  V  w 

w  v'w„  1 

l 

G 

G  / 

□ 


'  ( 

Q.  w 

+  ♦ 

*  M 

w„  Q  I 

l 

r 

G 

G 

G 

G  J 

r 


where 


a 

r 

D 

V4 

w ,  WQ 

♦.  »G 

M,MG 

0,Qq 


plate  domain 
plate  boundary 
flexural  rigidity 
biharmonic  operator 
deflections 
s 1  opes 

bending  moments  per  unit  length 
shear  forces  per  unit  length 


To  arrive  at  the  boundary  integral  equation,  wG  is  selected  as 
a  Green's  function. 
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3 . 2  Beam  Bending  Problem 

Butterfield  has  presented  a  reduced  form  of  Eq.  (3-1-1) 
for  the  one-dimensional  case  of  beam  bending,  which  is 
developed  in  the  discussion  that  follows.  As  a  first  step, 
Butterfield  multiplies  both  sides  of  the  governing  differential 
equation  for  a  beam  by  a  second  function  and  integrates  over 
the  beam’s  length,  giving 


El  V  w  w  dx 

Vj 


0 


A  w„  dx 

Cj 


(3-2-1 ) 


0 


where  A  is  the  distributed  load  per  unit  length,  I  is  the  area 
moment  of  inertia  of  the  beam  cross  section,  and  the  operator 
V4  is  defined  for  the  beam  case  as 


d 

dx 


Integrating  the  left  hand  side  of  Eq.  (3-2-1)  by  parts 
several  times,  Butterfield  develops  the  equation 


A  w. 


-  w  El  7 


‘WG  ) 


dx 


-WG  0  +  *G  M  "  MG  *  +  °G  w 


Recognizing  that  for  the  second  function 


E  1  V  WG  =  *G 

the  preceding  equation  can  be  rewritten  as 
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k  wG 

- 

W  k„ 

G 

dx 

= 

0 

-wQ  Q  +  *Q  M  -  MG  ♦  +  QG  w  (3-2-2) 

J0 


To  gain  some  physical  insight  into  Eq.  (3-2-2)  it  is 
useful  to  recall  Betti’s  reciprocal  work  theorem.  This  theorem 
states  that  for  an  elastic  structure  subjected  to  two 
independent  causes  (e.g.  loads),  the  total  work  done  by  the 
first  cause  in  moving  through  the  displacements  resulting  from 
the  second  cause  is  equal  to  the  total  work  done  by  the  second 
cause  in  moving  through  the  displacements  produced  by  the  first 
cause^^.  The  theorem,  as  it  applies  to  our  beam  bending 
problem,  can  be  better  understood  by  examining  the  examples  of 
Figure  3- 1 . 
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Figure  3  -  i  .  Examples  of  Betti’s  reciprocal  work  then  ■  i 

In  Figure  3-la,  the  product  of  load  and  displacement 
6xj  is  equal  to  the  product  of  load  Pj  and  displacement  6j^. 
Similarly,  in  Figure  3-lb,  the  product  of  moment  and  slope 
♦ij  is  equal  to  the  product  of  moment  Mj  and  slope  ♦ji- 

This  concept  can  be  applied  to  understand  Eq.  (3-2-2), 
with  the  i  cause  the  actual  loading  or  applied  moment  condition 
for  the  system  of  interest  and  the  j  cause  the  loading  or 
applied  moment  condition  of  the  second  function.  In  this 
equation,  the  A  term  represents  the  external  load  in  the  system 
of  interest,  which  when  multiplied  by  the  deflection  wq 
described  by  the  second  function,  gives  a  reciprocal  work  term. 
Similarly,  the  term  is  a  slope  for  the  second  function, 
which  when  multiplied  by  the  moment  M  for  the  system  of 
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interest  yields  another  work  term.  This  concept  applies  to  all 
the  terms  of  the  equation,  which  are  combined  according  to 
Betti's  theorm  of  reciprocal  work  to  yield  the  given  equality. 

For  Eq.  (3-2-2)  to  be  useful  for  solving  the  beam  bending 
problem,  Butterfield  chooses  as  the  second  function  a  Green’s 
function  which  is  based  on  a  singular  loading  condition;  that 
is  to  say,  a  point  load.  If  the  point  load  of  unit  magnitude 
is  applied  at  either  x  =  0  or  x  =  L,  the  second  intergral  term  on 
the  left  hand  side  of  Eq.  (3-2-2)  takes  on  the  value  of  the 
boundary  condition  of  the  beam  of  interest.  Usually,  some 
attention  should  be  given  to  cases  where  the  integrals  become 
s ingular. 

This  yields  two  equations.  However,  in  general,  beam 
bending  problems  involve  four  unknown  boundary  conditions.  For 
example,  if  simple  supports  are  specified  on  both  ends,  the 
moments  and  deflections  at  the  two  ends  are  known  to  be  zero, 
while  the  shears  and  slopes  are  not  known.  Therefore,  two 
additional  equations  are  required. 

To  obtain  two  more  equations,  derivatives  of  the  first  two 
equations  are  taken.  This  gives  an  additional  Green’s 
function.  Boundary  conditions  can  now  be  obtained  by  solving 
the  four  equations  simultaneously.  Knowing  the  boundary 
conditions,  the  beam  problem  can  be  solved  for  interior  points 
using  Eq.  (3-2-2)  and  appropriate  derivatives  with  the  point 


30 


load  positioned  at  the  location  of  interest. 


The  formulation  for  the  beam  bending  problem  outlined 
above  is  a  direct  approach  based  on  Green’s  identity. 
Butterfield  also  develops  a  boundary  integral  equation  using 
the  more  intuitive  indirect  method.  For  the  circular  plate 
formulation  that  follows  in  Chapter  4,  only  the  direct  method 
is  used . 
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CHAPTER  4 


DEVELOPMENT  OF  THE  FORMULATION  FOR  CIRCULAR  PLATES 


This  chapter  presents  a  general  method  for  solving 
circular  plate  bending  problems  in  the  elastic  range  using 
boundary  integrals.  The  resulting  integral  equations  are 
analagous  to  those  developed  by  Butterfield  for  the  one¬ 
dimensional  beam  problem,  which  are  discussed  in  Chapter  3. 
Treatment  of  non-linear  behavior  is  discussed  in  Chapter  5. 

4 .  I  Description  of  Plate  Geometry 

Figure  4  1  depicts  a  symmetrically  loaded  annular  plate, 

and  includes  notation  that  will  be  referred  to  throughout  the 
development  of  the  integral  equation  and  the  discussion  of  the 
selected  Green’s  functions. 

Important  symbols  are  defined  as  follows: 
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A 


.  2  Development  of  Integral  Equation 


Though  the  laterally  loaded  circular  plate  bending  problem 
is  two-dimensional,  symmetry  reduces  the  problem  to  essentially 
a  one-dimensional  type  problem.  Therefore,  it  is  expected  that 
the  integral  equation  developed  by  Butterfield  for  the  beam 
case  can  be  adapted  to  the  case  of  the  circular  plate.  This 
equation  (Eq.  (3-2-2))  is  repeated  below: 


I 

\  wG 

- 

w  \G 

dx 

= 

'o 

L 

~wa  Q  *  *g  M  ~  mg  *  +  qg  w 

Jo 

To  obtain  some  indication  of  what  the  circular  plate 
integral  equation  will  look  like,  a  non-rigorous  method  is 
first  used  to  adapt  Eq.  (3-2-2)  to  the  circular  plate  case, 
followed  by  a  more  rigorous  mathematical  development. 

In  the  non-rigorous  approach,  a  first  step  is  to  replace 
the  integral  on  the  left  hand  side  with  an  area  integral  and 
the  line  load  per  unit  length  \  with  the  load  per  unit  area  q. 
Next,  adjustments  are  made  to  reciprocal  work  terms  on  the 
boundaries  'right  hand  side  of  Eq.  (3-2-2)).  Specifically,  the 
shear  terras  are  replaced  with  the  product  of  the  radial  shear 
per  unit  length  and  the  total  length  (  2  *  r  )  ,  observing  that  for 
symmetrical  loading,  the  tangential  shear  Qg  is  0.  Moment 
terms  are  replaced  with  radial  moments,  since  there  are  no 
tangential  moments  on  the  boundaries.  These  moments  must  also 
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be  multiplied  by  the  length  (2ir)  to  give  total  moment.  The 
resulting  expression  is 


2  *  a 

0  b 


q  w. 


w  q_  r  dr  d8  = 

Vj 


2  *  r  j 

—  w  Q  + 

♦  _  M 

M.  ♦  +  Q_  w  I 

G  r 

G  r 

Gr  Gr  J 

which  reduces  to 


a 

(  q  wG  -  w  qQ  )  r  dr  = 
b 


r  I 

-w_  0  + 

M  - 

♦ 

Q„  w  ) 
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G  r 

G  r 

Gr 

Gr  j 

So  far,  it  has  only  been  asserted  that  the  integral 
equation  for  the  circular  plate  case  should  have  a  form  similar 
to  Eq.  (4-2-2).  To  obtain  the  exact  integral  equation,  a 
direct  method  is  used. 


As  a  starting  point  for  the  direct  method,  consider  the 
following  expression: 


U 


3  d  !w  d  *  w 

- _G  __  r  dr 

b  dr'  dr‘ 


This  expression  can  be  written  slightly  differently  as  V 


A 
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a  ,  2  ,2 

d  w  d  w„ 


dr  ‘  dr 


r  dr 


Integrating  both  expressions  by  parts  yields 
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d  2w 
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8 

1  d  2w 
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b  dr 
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dw 

d  2  WG  1 

a 
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^G 

dr 
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r  dr 2 

dr  3 

r  dr 


Recalling  Eq.  (2-3-2),  it  is  possible  to  rewrite  these 
expressions  in  terms  of  shear  Qr 


dw„  d  w 

Cl 

dr  dr  2 


dwG  Qr 
—  —  r  dr 

dr  D 


1  dw  dw 
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dw  d  w , 
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dw  QQr 
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1  dw  dw 
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r  dr  dr 


Setting  U  equal  to  V  and  canceling  some  terms  gives 
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dwG 

d  *  w 
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The  two  integral  terms  above  are  integrated  by  parts,  with 
the  operator  of  Eq .  (2-3-5)  used  to  simplify  the  expressions 
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which,  when  substituted  into  the  previous  equation,  yield 


dw  d  w 
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dr  dr 
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After  some  simplification,  this  expression  becomes 


7  w  r  w. 
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The  last  equation  bears  a  close  resemblance  to  Eq.  (4-2- 
2),  except  for  sign  differences  and  the  fact  that  terms 
involving  the  second  derivative  of  deflection  have  yet  to  be 
replaced  by  terms  involving  moment.  To  accomplish  the  latter, 
the  following  expression  is  added  to  the  third  terra  of  the 
right  hand  side  of  the  preceding  expression,  and  subtracted 
from  the  fourth  term 

v  dw  dw 

Ci 

r  dr  dr 


Through  Eq.  (2-1-7)  this  gives 


a 
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which  can  be  rewritten  using  Eq.  '2-3-4;  as 
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4-2  4  ' 


Remarkably,  with  the  exception 
attributable  to  the  difference 
expression  is  identical  to  Eq. 
based  on  the  integral  equation 


of  a  sign  difference  that  is 
in  sign  convention,  this 
(4-2-2),  which  was  asserted 
developed  by  Butterfield. 
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4 . 3  Selection  of  a  Green’s  Function 

As  is  the  case  for  the  beam  described  in  Chapter  3,  Eq. 

: 4-2-4)  is  useful  only  if  wg  is  chosen  as  a  Green’s  function 
with  certain  properties.  As  will  become  evident  later  in  this 
chapter,  the  Green’s  function  must  describe  a  ring  loaded  plate 
to  solve  circular  plate  bending  problems,  just  as  a  Green’s 
function  describing  a  point  loaded  beam  was  used  to  solve  the 
beam  bending  problem. 

To  solve  the  beam  problem,  Butterfield  choses  a  Green’s 
function  describing  an  infinitely  long  beam  with  a  point  load. 
It  is  not  essential  that  a  Green's  function  for  an  infinite 
beam  be  used.  In  fact,  as  one  would  expect  from  Betti’s 
reciprocal  work  theorem,  a  Green’s  function  describing  a  finite 
geometry  will  work  as  well,  provided  the  function  also 
describes  the  response  to  a  point  load. 

Applying  Butterfield’s  method  to  the  circular  plate  case, 
a  Green’s  function  representing  the  ring  loaded  plate  shown  in 
Figure  4-2  is  adopted.  The  Green’s  function  and  associated 
derivatives  are  included  in  Appendix  A.  It  should  be  noted 
that  the  outer  and  inner  radii  of  the  Green’s  function  <c  and 
d 1  need  not  coincide  with  the  outer  and  inner  radii  of  the 
actual  plate  < a  and  b , . 
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That  a  simple  support  is  chosen  for  the  outer  plate  radius  is 
not  important;  a  plate  clamped  at  both  the  inner  and  outer 
radii  would  also  have  been  suitable.  What  is  important  is  that 
the  loading  condition  for  the  Green’s  function  case  be  a  ring 

load.  To  simplify  calculations,  a  ring  load  of  magnitude  1  is 
used. 

4 . 4  Solution  of  Integral  Equation 

Having  selected  a  Green’s  function,  the  next  step  is  to 
obtain  a  form  of  Eq.  (4-2-4)  suitable  to  solve  the  actual  plate 
bending  problem.  Rearranging  terms,  Eq.  (4-2-4)  can  be 


r ewr i 1 1  en  as 
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The  ring  load  may  be  expressed  in  terms  of  the  dirac  delta 
function  S  (  r  ,  r  0  ) 


0  7  w. 


*(r,ro) 


Substituting  this  expression  into  Eq.  (4-4-1)  and  making  use  of 
the  following  property  of  a  dirac  delta  function 
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By  choosing  r0  at  radius  a,  and  then  at  radius  b,  the 
deflection  w  in  the  left  hand  side  of  Eq.  ,4-4-2)  is  the 
deflection  at  the  boundary,  which  is  either  a  known  or  unknown 
boundary  condition.  Hence,  two  equations  are  available  to 
solve  for  the  four  unknown  boundary  conditions. 
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To  obtain  two  more  equations,  Eq.  (4-4-2)  is 


differentiated  with  respect  to  r0.  This  yields  the  second 
Green’s  function,  which  is  evaluated  with  rQ  at  a  and  then 
again  at  b.  The  second  Green’s  function  corresponds  to  the 
slope  of  the  plate,  which  at  r=a  and  r=b  is  again  either  a 
known  or  an  unknown  boundary  condition.  There  is  thus  now  a 
complete  set  of  four  equations  to  solve  for  the  four  unknown 
boundary  conditions. 


Eq.  (4-4-2)  can  be  rewritten  as 
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Differentiating  this  expression  with  respect  to  rD  gives 
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Eqs.  (4-4-3)  and  (4-4-4),  evaluated  both  for  rQ=a  and  r0=b, 
give  rise  to  a  system  of  four  equations  that  when  solved 
simultaneously,  yield  the  four  unknown  boundary  conditions.  In 
matrix  form,  these  equations  can  be  written  as 
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Using  Gauss  elimination,  the  preceding  system  of  equations  is 
solved  for  the  four  unknown  boundary  conditions. 


The  final  step  is  to  calculate  deflection,  slope,  moment 
and  shear  across  the  plate.  Deflection  and  slope  are  found 
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using  Eqs.  (4-4-3)  and  (4-4-4)  by  substituting  in  the 
calculated  boundary  conditions  and  evaluating  the  expressions 
at  the  value  of  rQ  of  interest.  Moment  and  shear  are 
determined  using  Eqs.  (2-1-7)  and  (2-3-2),  which  require 
calculation  of  the  second  and  third  derivatives  or  w  with 
respect  to  rQ.  Differentiating  F. q.  (4-4-4)  results  in  the 
desired  expressions 
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making  a  complete  solution  for  the  circular  plate  bending 
problem  in  the  elastic  range  possible. 
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CHAPTER  5 


NON-LINEAR  BEHAVIOR 


This  chapter  develops  a  boundary  integral  formulation  for 
circular  plate  bending  problems  that  includes  non-linear 
material  behavior.  Specifically,  the  boundary  integral 
equation  developed  in  Chapter  4  for  the  elastic  plate  is 
modified  to  include  plastic  moment  terms  to  account  for  the 
plastic  behavior.  An  incremental  load  method  similar  to  that 
used  by  Moshaiov  and  Vorus  can  then  be  applied  to  arrive  at  a 
solution  to  the  plate  bending  problem  once  yielding  has  begun. 

5 .  1  Plastic  S t r ess -S t r a  i  n  Relationships 

When  considering  plasticity,  strain  can  be  thought  of  as 
having  two  components:  an  elastic  component  and  a  plastic 
component.  In  cylindrical  coordinates,  which  is  applicable  to 
the  two-dimensional  circular  plate  case,  the  total  strain  t  can 
be  expressed  in  terms  of  these  two  components  as 
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where  the  superscripts  e  and  p  refer  to  the  elastic  and  plastic 
strain  components,  respectively. 


Eq.  (2-1-4),  which  provides  the  relationship  between 
stress  and  strain  for  a  linearly  elastic  material,  can  thus  be 
rewritten  in  terms  of  the  total  strain  and  plastic  strain  as 
f  o 1  lows : 
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5 . 2  Plastic  Moment-Curvature  Relationships 

To  develop  the  moment-curvature  relationships  with 
plasticity,  Eq.  (5-1-2)  is  substituted  into  Eq.  (2-1-3)  to  give 
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Recalling  the  relationships  for  total  strain  of  Eq.  (2-1-6)  and 
performing  the  integration  over  the  plate  thickness  yield  the 
moment-curvature  relationship  for  a  circular  plate  with 
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plast  i  c  i  t  y 
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(5-2-3) 


5 . 3  Governing  Differential  Equation  With  Plasticity 

Chapter  2  develops  the  governing  differential  equation  for 
a  circular  plate  in  the  elastic  range  by  substituting  the 
elastic  raoinen t - curve t ure  relationships  into  the  elastic 
equilibrium  equation.  In  this  section,  a  similar  approach  is 
used  to  develop  the  governing  differential  equation  with 
plasticity.  The  difference  here  is  that  the  moment-curvature 
relationships  (Eq.  (5-2-3))  include  plastic  terms. 


Substituting  Eq.  (5-2-3)  into  the  equilibrium 
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(5-3-1  ) 


where  the  tilda  is  used  to  distinguish  the  fact  that  the  shear 
includes  plastic  effects.  Differentiating  this  expression  once 
with  respect  to  r  gives  the  fourth  order  governing  differential 
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equation  for  the  plate  with  plasticity 
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which  can  be  rewritten  using  the  operator  defined  in  Eq.  ; 2  —  3  - 
5  )  as 
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This  equation  is  identical  to  the  equilibrium  equation  for  the 
elastic  case  with  the  exception  of  the  terms  that  include 
derivatives  of  the  plastic  moments.  Note  that  these  terms  may 
be  thought  of  as  additional  pseudo  loads  that  must  be  combined 
with  the  actual  load  to  account  for  the  plastic  behavior. 


5 . 4  Boundary  Integial  Formulation  With  PI  as t  i  city 

In  the  preceding  chapter,  the  governing  differential 
equation  is  used  to  develop  a  boundary  integral  formulation  for 
a  circular  plate  in  the  elastic  case.  This  section  presents  a 
similar  approach  to  develop  a  boundary  integral  formulation 
for  the  plastic  case. 


As  a  starting  point,  Eqs.  (2-2-4)  and  (5-3-2)  are 
substituted  into  Eq.  (4-2-3)  to  give 
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This  equation  is  not  yet  in  a  form  that  is  usable  for 
solving  the  circular  plate  bending  problem.  In  the  first 
place,  there  are  several  terms  that  include  derivatives  of  the 
radial  or  tangential  plastic  moments,  which  are  not  generally 
known  across  the  plate.  Secondly,  the  shear  and  moment  terms 
Qr  and  Mr  evaluated  at  the  boundaries  are  not  the  actual  shears 
and  moments,  since  they  were  developed  based  on  elastic 
considerations  only. 


To  resolve  these  concerns,  the  integral  term  involving 
plastic  moments  in  Eq.  (5-4-1)  is  integrated  by  parts  as 
foil ows : 
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integrating  by  parts  a  second  time  and  collecting  terms  gives 


Using  the  above  relationship  and  choosing  a  Green’s  function 
that  describes  a  ring  loaded  plate  as  is  done  for  Eq.  (4-4-2) 
Eq.  (5-4-1)  becomes 
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Recall  that  Qr  and  Mr  represent  shear  and  moment  for  the 
elastic  case.  However,  they  are  not  the  actual  shear  and 
moment  for  the  plastic  case.  With  plasticity,  the  actual 
shears  and  moments  are  combinations  of  elastic  and  plastic 
terms.  From  Eq.  (5-3-1),  the  shear  relationship  with 
plasticity  is 


0 


r 


and  from  Eq.  (5-2-3),  the  moment  relationship  with  plasticity 
i  s 


=  M  +  M  p 
r  r  r 

where,  as  with  the  case  for  shear,  the  tilda  is  used  to 
distinguish  the  fact  that  the  moment  terra  includes  plastic 
effects.  When  substituted  into  the  preceding  expression,  these 
expressions  give  the  desired  boundary  integral  formulation  for 
*he  circular  plate  bending  problem 
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By  evaluating  Eq.  (5-4-2)  with  the  ring  load  positioned 
both  at  r0-a  and  r0=b,  two  equations  are  obtained  to  solve  for 
the  unknown  boundary  conditions.  Eq.  (5-4-2)  is  differentiated 
with  respect  to  r0  to  obtain  a  second  Green’s  function.  This, 
via  Eq.  (2-1-1),  gives  an  expression  for  slope  of  the  actual 
plate.  The  second  Green's  function  is  evaluated  with  the  ring 
load  positioned  at  r0=a  and  rQ-b,  to  provide  the  remaining  two 
equat  ions . 


Differentiating  Eq.  (5-4-2)  with  respect  to  r0  gives 
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In  matrix  form,  the  system  of  equations  to  be  solved  for 
the  four  unknown  boundary  conditions  for  the  non-linear  case 
can  be  written  using  the  notation  of  Eq.  (4-4-5)  as  follows: 
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Observe  that  Eq.  (5-4-4)  is  identical  to  Eq.  (4-4-5) 
except  for  the  additional  moment  terms  in  the  last  matrix  of 
Eq.  (5-4-4).  These  additional  terms  account  for  the  plastic 
behav ior . 


To  find  moments  and  shears  across  the  plate,  expressions 
for  the  second  and  third  derivatives  of  deflection  with  respect 
to  r0  are  required.  These  expressions  follow: 
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CHAPTER  6 


NUMERICAL  SOLUTION 


This  chapter  describes  the  numerical  procedure  that  was 
Used  to  solve  the  circular  plate  bending  problem  in  both  the 
elastic  and  plastic  ranges.  The  incremental  load  approach  used 
in  the  plastic  range  is  first  described.  A  brief  description 
of  the  computer  program  that  was  developed  based  on  the 
equations  of  Chapters  4  and  5  follows,  including  a  discussion 
of  some  aspects  of  numerical  methods  used  in  this  program. 

6  •  1  Incremental  Load  Method 

Once  plastic  deformation  starts  to  occur  in  the  plate,  the 
linear  relation  between  stress  and  strain  is  not  valid.  As 
discussed  by  Mendelson*'®'*,  strains  in  the  plastic  range  are  no 
longer  uniquely  determined  by  the  stresses,  and  in  fact  depend 
on  the  loading  history.  Therefore,  it  is  not  possible  to  use 
the  relationships  developed  in  Chapter  5  to  directly  arrive  at 
a  plate  bending  solution  given  a  loading  condition  outside  the 
elastic  range.  Rather,  an  incremental  approach  must  be  used 
wherein  the  load  is  increased  in  small  steps  once  yielding 
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occurs  at  any  point  in  the  plate,  and  the  complete  stress  state 
for  the  plate  is  determined  before  load  again  is  increased. 

The  incremental  load  method  used  for  this  analysis  is 
based  on  that  adopted  by  Moshaiov  and  Vorus,  and  is  frequently 
seen  in  finite  element  solutions  of  structural  analysis 
problems.  The  principal  steps  of  this  method  are  outlined 
below: 

1.  Using  the  plate  parameters  and  loading  conditions  for 
the  bending  problem  at  hand,  an  elastic  solution  is  obtained 
with  the  equations  developed  in  Chapter  4. 

2.  The  load  at  which  the  onset  of  yielding  will  occur  in 
the  plate  is  calculated  based  on  the  stress  state  determined  in 
Step  1.  The  selected  yield  criterion  is  the  von  Mises 
criterion,  which  for  the  symmetrically  loaded  circular  plate 
case  has  the  form: 


where  oe  is  the  equivalent,  or  effective,  stress  which 
represents  the  von  Mises  yield  surface.  Yielding  will  occur 
when  oe  equals  or  exceeds  the  uniaxial  yield  stress  { S  y )  . 

3.  Unknown  boundary  conditions  at  yielding  are 
determined  using  Eq.  (5-4-4),  with  plastic  moments  initially 
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set  equal  to  zero.  Total  strains  at  predetermined  integration 
points  are  calculated  using  Eqs .  (2-1-6),  (5-4-3)  and  (5-4-5). 

From  Eq.  (2-1-4),  the  stress  state  at  the  integration  points  at 
yielding  is  determined,  and  is  stored  to  be  later  updated  as 
load  is  increased.  Also  stored  for  later  updating  are  the 
deflections,  slopes,  moments  and  shears  across  the  plate. 

4.  The  load  which  caused  yield  is  then  increased  by  a 
small  incremental  amount.  The  increments  of  the  total  strains 
resulting  from  this  incremental  load  are  calculated  using  Eq. 
(2-1-6)  in  an  incremental  form,  from  which  elastic  strains  can 
be  calculated  by  means  of  Eq.  (5-1-1). 

5.  Using  Eq.  (2-1-4),  the  elastic  stress  increment 
corresponding  to  the  applied  incremental  load  is  calculated. 
This  incremental  stress  is  added  to  the  stress  stored  in  Step  3 
to  give  the  total  stress  at  the  end  of  the  load  increment.  A 
check  is  made  throughout  the  plate  using  the  yield  criterion  of 
Step  2  to  determine  where  yielding  has  occurred. 

6.  In  those  plate  regions  where  yielding  has  occurred, 
the  plastic  strain  increment  (that  is,  the  plastic  strain  due 
to  the  incremental  load)  is  calculated.  Only  materials  that 
exhibit  e 1  as t ic-per f ec 1 1 y  plastic  behavior  as  depicted  in 
Figure  6-1  have  been  examined.  However,  the  method  can  be  used 
with  other  material  behavior,  such  as  strain  hardening.  For 
materials  exhibiting  e 1  as t i c-per f ect 1 y  plastic  behavior,  the 
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plastic  strain  increment  Sc p  is  the  total  incremental  strain 
determined  from  Eq.  (2-1-6).  This  is  shown  in  Figure  6-1. 

7.  The  plastic  moments  are  determined  using  Eq.  (5-2-2). 

8.  Steps  3  through  6  are  repeated,  with  the  plastic 
moments  calculated  in  Step  7  used  as  an  input  to  Eq.  (5-4-4). 
Step  7  is  repeated  and  the  plastic  moments  compared  to  the 
moments  calculated  previously  in  Step  7.  Iterations  are 
performed  as  necessary  until  these  values  converge,  according 
to  a  predetermined  convergence  criteria. 

9.  After  convergence  of  plastic  moments  is  achieved  in 
Step  8,  stresses,  plastic  strains,  deflections,  slopes,  moments 
and  shears  throughout  the  plate  are  updated. 

10.  Another  increment  of  load  is  applied,  and  Steps  3 
through  9  above  are  repeated,  using  the  plastic  moment  from  the 
previous  load  step  as  an  initial  estimate  for  plastic  moment. 
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6  •  2  Computer  Program  D escript ion 

Computer  programs  written  in  the  Fortran  77  language  have 
been  developed  to  solve  the  circular  plnte  bending  problem  and 
are  included  as  Appendix  B.  Three  programs,  supported  by  nine 
subroutines,  are  used. 

The  programs  INPUT  and  LOAD  create  data  files  which  are 
used  by  the  program  MAIN  to  solve  the  bending  problem.  INPUT 
and  LOAD  prompt  the  user  for  information  describing  the  problem 
to  be  solved,  such  as  material  properties,  plate  geometry, 
boundary  conditions  and  plate  loading  conditions.  In  addition, 
certain  adjustable  parameters  are  input,  which  include  step 
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sizes  for  the  load  increments,  number  of  integration 
increments,  and  the  desired  percentage  for  convergence  of  the 
plastic  moment  increments.  The  program  MAIN  determines  the 
plate  bending  solution,  both  in  the  elastic  and  plastic  ranges. 

A  flow  diagram  outlining  the  major  elements  and  overall 
logic  path  of  this  program  is  shown  in  Figure  6-2. 
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Figure  6-2.  Computer  program  flow  diagram. 
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6 . 3  Numerical  Methods 


Eqs.  (5-2-1),  (5-4-4),  (5-4-5)  and  (5-4-6)  require 
integrations  to  be  performed  across  the  plate.  These 
integrations  are  accomplished  by  a  trapozoidal  rule  numerical 
integration  scheme.  It  was  found  that  10  increments  across 
both  the  plate  radius  and  the  plate  half  thickness  gave 
satisfactory  results. 

To  ensure  convergence  of  plastic  moments,  a  root  mean 
square  average  of  the  the  plastic  moments  at  each  station 
across  the  plate  is  calculated  and  compared  to  the  root  mean 
square  average  from  the  previous  iteration.  An  agreement  of 
less  than  0.1%  was  used  for  the  analysis  work  for  which  results 
appear  in  Chapter  7. 
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CHAPTER  7 


RESULTS  AND  DISCUSSION 


This  chapter  presents  results  of  analyses  performed  on 
three  different  plate  configurations  using  the  computer  program 
described  in  the  preceding  chapter.  Results  in  the  elastic 
range  for  all  three  cases  are  compared  with  the  analytic 
solution  from  Roark<^>,  and  for  one  case  in  the  plastic:  range 
where  published  results  using  another  solution  method  are 
available.  Descriptions  of  the  case  studies  are  given  in 
Sections  7.1  through  7.3.  A  discussion  of  the  results  is  given 
in  Section  7.4. 

7 . 1  Simply  Supported  Elasto-Plastic  Annular  Plate 

The  annular  plate  shown  in  Figure  7-1  is  first  examined, 
The  plate  is  simply  supported  at  both  the  inner  and  outer 
radii,  and  is  subjected  to  a  uniform  lateral  load.  A  material 
exhibiting  elastic-perfectly  plastic  behavior  as  depicted  in 
Figure  6-1  is  used.  The  yield  stress  (Sy)  is  16  ksi  and  the 
Young’s  modulus  is  10xl03  ksi.  Material  properties  and  plate 


dimensions  were  selected  to  allow  comparison  with  results 


obtained  using  a  finite  element  method  by  Armen  et  a  1 ^  ^  _ 

Plate  deflection  at  yielding  is  shown  in  Figure  7-2,  which 
also  includes  results  from  the  analytic  solution  given  by  Roark 
for  comparison.  Plate  deflections  for  the  present  solution  at 
several  loads  in  the  plastic  range  are  shown  in  Figure  7-3,  as 
well  as  results  obtained  by  Armen  et  al  using  a  finite  element 
method  of  analysis. 

Figure  7-4  shows  the  plastic  zones  in  the  plate  at  the 
three  load  conditions  depicted  in  Figure  7-3.  Results  obtained 
by  Armen  et  al  are  also  included. 


Sy  =  16  KSI 
E  =  10X103  KSI 

V  =  0.24  1 

i  10,0  IN 


Figure  7  -  I  .  Simply  support*?'!  annular  r  •  a  t  p  r  o  L>  i  m 
de  s cription. 
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7 . 2  Clamped  Elasto-Plastic  Annular  Plate 

The  annular  plate  shown  in  Figure  7-5  is  next  examined. 

The  plate  is  clamped  at  both  the  inner  and  outer  radii,  and  is 
again  subjected  to  a  uniform  lateral  load.  A  material 
exhibiting  the  same  properties  as  described  in  Section  7.  I  is 
used.  For  the  clamped  annular  plate  case,  results  obtained 
using  other  methods  could  not  be  found  for  the  plastic  range. 
Consequently,  material  properties  and  plate  dimensions  were 
selected  to  allow  a  qualitative  comparison  with  results  for  the 
simply  supported  plate  of  Section  7.1 

Plate  deflections  at  the  onset  of  yielding  are  shown  in 
Figure  7  -  6 ,  and  at  several  loads  in  the  plastic  range  in  Figure 
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7-7  . 


The  deflection  at  yielding  using  the  analytic  solution 


given  by  Roark  is  shown  in  Figure  7-6  for  comparison.  Finally 
Figure  7-8  shows  the  plastic  zones  in  the  plate  at  the  three 
load  conditions  depicted  in  Figure  7-7. 
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figure  7-8.  Plastic  zones  in  clamped  annular  piste. 

7 . 3  Simply  Supported/Guided  E 1  as t o-P 1  as t  i  c  Annular  Plate 

The  third  annular  plate  to  be  examined  is  shown  in  Figure 
7-8.  The  plate  is  simply  supported  at  the  outer  radius,  with 
the  inner  radius  guided.  The  plate  is  again  subject  to  a 
uniform  lateral  load.  A  material  exhibiting  the  same 
properties  described  in  Section  7-1  is  assumed. 

For  this  case,  results  obtained  using  other  methods  could 
again  not  be  found  for  the  plastic  range.  However,  a  solution 
for  a  continuous  simply  supported  circular  plate  has  been 
obtained  by  Moshaiov  and  Vorus,  and  offers  a  qualitative 


comparison  with  the  results  for  the  simply  supported/guided 


annu  lar  plate. 


Plate  deflection  at  the  onset  of  yielding  is  shown  in 
figure  7-10  and  at  several  loads  in  the  plastic  range  in  Figure 
7-11  for  the  present  analysis  method.  The  deflection  at 
yielding  obtained  using  the  analytic  solution  from  Roark  is 
shown  in  Figure  7-10  for  comparison.  Figure  7-12  shows  the 
plastic  zones  in  the  plate  at  the  three  load  conditions 
depicted  in  Figure  7-11  for  the  present  solution.  Corres¬ 
ponding  results  from  Moshaiov  and  Vorus  for  the  continuous 
;■!  ate  of  Figure  7-13  are  shown  in  Figures  7-14  and  7-15. 
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7.4  D iscuss 1  on  of  Results 


The  results  presented  in  Sections  7.1  through  7.3  for  the 
elastic  range  are  uniformly  in  excellent  agreement  with  the 
results  using  the  analytic  solutions  given  by  Roark.  This  is 
as  should  be  expected,  since  the  elastic  range  solution  using 
the  present  method  is  based  on  a  closed  form  analytic  solution. 

Results  in  the  plastic  range,  while  reasonable,  are  not  in 
as  close  agreement  with  other  published  results.  The  only  case 
where  data  for  a  direct  comparison  could  be  found  is  the 
configuration  discussed  in  Section  7.1.  From  Figure  7-3,  the 
deflections  at  loads  of  652  psi  and  852  psi  for  the  present 
solution  are  in  excellent  agreement  with  the  finite  element 
results  of  Armen  et  ai. 

Results  for  the  annular  plate  of  Section  7.1  at  the  1052 
psi  load  are  not  as  encouraging.  Deflections  and  the  extent  of 
the  plastic  zone  are  greater  for  the  solution  obtained  by  Armen 
than  for  the  present  solution.  To  explain  this  difference,  the 
size  of  the  integration  and  load  increments  has  been  varied, 
but  it  shows  little  effect  on  the  results.  The  reason  behind 
this  difference  is  not  known.  However,  it  is  believed  that  it 
might  be  due  to  a  programming  error  either  in  this  report  or  in 
the  reference.  Also,  differences  could  arise  by  not  taking 
into  account  additional  integration  terms  arising  from  the 
singularity  at  r=rQ,  as  discussed  in  Appendix  A. 
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Finally,  it  is  of  interest  to  compare  the  results  for  the 
simply  supported/guided  plate  of  Figure  7-9  to  the  results  for 
a  s imp l y-s uppor t ed  circular  plate  obtained  by  Moshaiov  and 
Vorus.  It  is  recognized  that  the  two  cases  are  not  equivalent. 
However,  one  would  expect  the  two  plates  to  behave  globally  in 
a  similar  fashion.  This  is  in  fact  the  case,  with  both 
deflections  and  plastic  zones  showing  reasonably  good  agreement 
in  view  of  the  differences  between  the  two  configurations. 


CHAPTER  8 


SUMMARY  AND  CONCLUSIONS 


8 . 1  Summary 

This  thesis  has  developed  a  formulation  for  solving 
ax  i  symme t r i ca 1 1 y  loaded  annular  plate  bending  problems  using 
boundary  integrals.  This  particular  formulation  is  unique  in 
that  it  treats  the  annular  plate  as  a  one-dimensional  problem, 
using  "ring"  type  Green’s  functions  to  determine  unknown 
boundary  conditions  and  arrive  at  a  plate  bending  solution. 

The  formulation  allows  a  numerical  incremental  load  method  to 
be  used  in  the  plastic  range  for  the  treatment  of  non-linear 
behavior.  Results  in  the  elastic  range  show  excellent 
agreement  with  results  obtained  using  a  conventional  analytic 
solution.  In  the  plastic  range,  results  are  in  reasonable 
agreement  with  those  obtained  using  a  finite  element  method. 

8 . 2  Cone  1  us i ons 

This  thesis  treats  the  two-dimensional  problem  of  analysis 
of  ax i s yrame t r i ca 1 1 y  loaded  annular  plates  as  a  one-d iraens i ona  1 
problem.  Using  a  method  of  solving  one-dimensional  problems 
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with  boundary  integrals  developed  by  Butterfield,  the  thesis 
demonstrates  that  a  closed  form  solution  for  the  annular  plate 
problem  can  be  obtained.  It  further  demonstrates  that  this 
approach  is  suitable  for  the  use  of  an  incremental  load  method 
to  obtain  a  solution  in  the  plastic  range. 

The  formulation  developed  in  this  thesis  is  advantageous 
over  a  conventional  analytic  solution  in  that  it  allows  a  wide 
range  of  problems  with  different  loading  and  boundary 
conditions  to  be  solved  using  a  single  algorithm.  Further,  the 
simplicity  of  the  approach  is  useful  from  an  educational 
standpoint  in  the  teaching  of  boundary  element  methods. 

8.  .'J  Recommendations  for  Future  Work 

There  remains  room  for  much  additional  work  in  this  area. 
Some  suggestions  follow: 

1.  A  formulation  for  ax i symmet r i ca 1 1 y  loaded  continuous 
circular  plates  remains  to  be  developed.  The  work  presented 
here  lays  a  foundation  for  the  continuous  plate  formulation. 
However,  extending  this  method  to  the  continuous  plate  case  may 
require  some  modifications. 

2.  In  developing  a  formulation  for  the  simple  beam 
problem  using  boundary  integrals,  Butterfield  uses  both  a 
direct  and  an  indirect  method,  and  shows  that  they  yield  the 
same  result.  Only  the  direct  method  is  presented  in  this 
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thesis.  The  indirect  method  is  a  more  intuitive  approach  than 
the  direct  method,  and  offers  advantages  in  furthering  the 
understanding  of  a  boundary  integral  solution.  It  would, 
therefore,  be  worthwhile  to  pursue  an  indirect  method  as  well 
as  a  direct  method  in  developing  a  boundary  integral 
formulation  for  circular  plate  geometries. 

3.  The  Green’s  functions  selected  for  this  anlysis 
involve  many  terms  and  are  awkward  from  a  computational 
standpoint.  More  thought  needs  to  be  given  to  selecting  a 
"ring"  type  Green’s  function  that  is  simpler  and  therefore 
offers  advantages  in  simplifying  calculations  and  reducing 
computational  time. 

4.  For  the  plastic  range,  the  lack  of  close  agreement 
between  the  results  presented  in  this  work  and  those  obtained 
by  Armen  et  al  using  a  finite  eleroet  method  need  to  be  better 
understood.  Additional  comparative  data  should  be  found  or 
developed  to  increase  confidence  in  the  results  presented  in 
this  work. 
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APPENDIX  A 


SELECTED  GREEN’S  FUNCTIONS 


This  Appendix  presents  the  selected  Green’s  functions  and 
required  derivatives  that  are  used  to  solve  the  plate  bending 
problems  in  this  analysis.  The  functions  which  mathematically 
describe  the  deflection,  slope,  radial  and  tangential  moments, 
and  shear  are  obtained  from  Roark. 

A . 1  Description  of  Geometry  and  Sign  Conventions 

As  discussed  in  Chapter  4,  the  first  Green’s  function 
selected  for  this  analysis  describes  the  response  of  annular 
plate  that  is  simply  supported  at  the  outer  radius  and 
unsupported  along  the  inner  radius  to  a  ring  load.  Figure  A-l 
depicts  this  plate  configuration  and  important  notation. 
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Figure  A-l.  Representation  of  Green's  function  used  in 
the  analysis. 

Roark  uses  a  sign  convention  with  deflection  w q  positive 
upward,  slope  positive  when  the  deflection  wq  increases 
positively  as  r  increases,  moment  Mr  positive  when  creating 
compression  on  the  top  surface  of  the  plate,  and  the  shear 
force  Qcr  positive  when  acting  upward  on  the  inner  edge  of  an 
annular  section.  Subscripts  c  and  d  refer  to  the  radial 
pos i t  1  on  . 

This  sign  convention  differs  from  the  sign  convention  of 
Timoshenko,  which  has  been  adopted  for  this  analysis. 
Accordingly,  adjustments  must  be  made.  Specifically, 
Timoshenko  takes  deflection  to  be  positive  downward  and  shear 
to  be  positive  when  acting  downward  on  the  inner  edge  of  the 


plate,  such  that  the  signs  for  these  items  as  given  by  Roark 
mus  t  be  changed . 

For  slope,  Roark  uses  a  convention  opposite  to  that  of  Eq. 
(2-2-1)  and  the  convention  used  by  Timoshenko.  This 
difference,  combined  with  the  difference  in  defining  the  sign 
of  deflection,  effectively  cancel  each  other,  so  the  sign  of 
slope  from  Roark  does  not  change  for  this  analysis.  Finally, 
for  moments,  Roark  and  Timoshenko  use  the  same  sign  convention, 
so  no  sign  changes  for  moment  terms  are  necessary. 

The  formulas  in  the  sections  which  follow  include  a  number 
of  general  plate  functions  and  constants  that  are  not  depicted 
in  Figure  A  - 1 .  Functions  L  and  G  and  their  derivatives  are 
included  in  Section  A. 7;  constants  F  and  C  are  included  in 
Section  A. 8. 

A . 2  Singulari ties 

It  is  not  possible  to  evaluate  the  Green’s  function  and 
associated  derivatives  at  the  exact  radial  location  where  the 
ring  load  is  applied,  due  to  singularities.  To  avoid  this 
problem  on  the  boundaries,  the  ring  load  is  positioned  a  small 
distance  inside  the  outer  plate  radius  and  outside  the  inner 
plate  radius  (0.00001  inches  for  the  results  presented  in 
Chapter  7).  During  final  review  of  this  thesis,  it  was  noted 
’hat  because  of  the  singularity  at  r=rov  the  domain  integrand 
involving  plastic  moments  of  Eqs.  (5-4-2)  through  (5-4-6) 


includes  singularities.  Therefore,  the  limiting  values  should 
be  explored  and  a  correction  applied  in  way  of  the  singular¬ 
ities.  Insufficient  time  was  available  to  repeat  the  analysis 
with  the  applied  correction  to  see  its  effect  on  the  results. 
Note  that  this  correction  affects  only  the  solution  in  the 
plastic  range . 


A . 3  Deflection.  Slope.  Moment  and  Shear 

The  following  expressions  describe  deflection  wq ,  slope 
♦  q  ,  radial  moment  Mqj.,  and  shear  Qq  corresponding  to  the  first 
Green ’ s  f  unct ion : 
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D  =  plate  flexural  rigidity  (Eq.  (2-1-8 
P_  =  magnitude  of  ring  load 

Ij 

0 

<r-rQ>  =  singularity  function 


The  expression  <r-ro>0  represents  a  singularity  function. 
The  function  is  equal  to  zero  if  r<rQ.  If  r>rQ,  the  brackets 
become  like  any  other  brackets.  Hence,  for  r>rD,  <  r  -  r  0  >  0  is 
(r-r0)  to  the  power  of  0,  which  equals  to  1. 


A . 4  First  Derivative  of  Deflection.  Slope,  Moment  and  Shear 
The  expressions  which  follow  describe  the  first  deriv¬ 
atives  with  respect  to  rQ  of  wq,  $q,  Mqj-  and  Ocr-  The  first 
derivative  of  wq  corresponds  to  the  second  Green’s  function 
required  by  the  analysis. 
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A • 6  Third  Derivative  of  Deflection.  Slope.  Moment  and  Shear 
The  following  expressions  describe  the  third  derivatives 
with  respect  to  rQ  of  wG  ,  *G,  MGr  and  QGr 
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A . 8  Plate  Constants  f 


F 


F 


C 


and  C 
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function  l  and  MQr)  are  available  and  have  been  presented  in 
Section  A. 3.  It  is  further  noted  that  Eqs.  (5-4-2)  through  15- 
4-6)  require  derivatives  of  the  above  two  expressions  with 
respect  to  r0.  This  is  accomplished  via  the  expressions  for 
derivatives  of  and  M(jr  with  respect  to  rQ  included  in 
Sect  ons  A. 5  through  A. 7  above. 
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APPENDIX  B 


COMPUTER  PROGRAMS 


This  Appendix  contains  the  computer  programs  developed  to 
perform  the  analysis  work.  A  flow  chart  showing  the  overall 
program  structure  is  presented  ir.  Chapter  6. 

The  Appendix  is  organized  with  the  programs  INPUT  and  LOAD 
presented  first,  followed  by  the  program  MAIN.  The  supporting 
subroutines  for  the  program  MAIN  follow  that  program. 


PROGRAM  INPUT 

THIS  PROGRAM  ALLOWS  THE  USER  TO  DEFINE  THE  PROBLEM  OF  INTEREST 
AND  TO  SET  CERTAIN  ADJUSTABLE  PARAMETERS  SUCH  AS  NUMBER  OF  DEPTH 
INTEGRATION  INCREMENTS.  LOADING  CONDITIONS  ARE  INPUT  USING  THE 
PROGRAM  LOAD.  TO  RUN,  THIS  PROGRAM  MUST  ACCESS  THE  FILES  GEOM.DA 
BC.DAT.BCPOS.DAT  AND  ADJ . DAT 


PROGRAM  INPUT 

IMPLICIT  RE AL*8  iA-H.O-Z) 

CHARACTER* 1  B FLAG  1 1 4 )  , B FLAG2 ( 4 ; ,  RESP 
CHARACTER*19  LB  L ( 8 ) , LG EOM ( 9  > 
CHARACTER*25  L AD J ( 6 ) 

HE AL*8  ADJ(7) , BC(4) .GEOMC9) ,NU 

INTEGBH  POS( 8) , BCPOSf 8  '  ,M, N 


INPUT  MATERIAL  CONSTANTS  AND  PLATE  GEOMETRY.  THE  FLEXURAL 
R IG ID  T  T Y  IS  CALCULATED  BASED  ON  THE  INPUT  PARAMETERS  N  IF  THE 
TOTAL  NUMBER  OF  PARAMETERS 


N  -  9 

LGEOMM  = ’OUTER  RADIUS  A  =  ’ 

LGBOM; 2 ) = ’ INNER  RADIUS  B  =  ’ 

LGBOM{ 3 )=’ POISSON  RATIO  NU  =  ' 

LGBOM( 4 ) = ’ PLATE  CONSTANT  D  =  ’ 

LGEOM(5 ) =  ’ THICKNESS  T  =  * 

LGBOMC6) *' YIELD  STRESS  =  ’ 

LGK0M'7\= ’YIELD  STRAIN  -  ’ 

LGEOM, 8' = ’ ULTIMATE  STRESS  =  * 

LG EOM 1 9  - 'ULTIMATE  STRAIN  =  * 

WRITE’  *  .  1  ) 

1  FORMAT  '  ,  6X  ,  ’  WOULD  YOU  LIKE  TO  SE£  THE  CURRENT  PLATE  DIMES'-  ! 
*  ,  /  ,  2  X  ,  ’AND  MATERIAL  CONSTANTS  Y  N  °  * 

R E  A D  f  * . 30  RESP 

IF  RESP  . EO.  ' Y ’  THEN 
WRITE  *,2 


95 


2  FORMAT! // ,  2X ,  ’ THE  CURRENT  VALUES  FOLLOW.  NOTB  THAT  A  ZBBO', 

♦  /,2X, ‘VALUE  IS  SHOWN  FOB  D  WB8N  T  AND  B  ARB  GIVEN  AND  VICE’. 

♦  / , 2  X ,  ’VERSA 

OPEN { R ,  F I  LB  = ’ G  BOM . D  AT  *  , STATUS = ’OLD ' ) 

WRITER*, 30)  *  • 

DO  8  1=2, N 

R8AD(6,70)  GBOM(I) 

WRITE ( * , 6 )  LG BOM ( I ) , GBOM ( I ) 

6  FORMAT(2X. A22, F20. 10) 

8  CONTINUE 
ClOSB(6) 

END  I  F 
C 

WRITE!*, 9) 

9  FORMAT! /', 6X DO  YOU  WISH  TO  INPUT  OR  CHANGE  GEOMETRY  OR’, 

♦  /,2X, ’MATBRIAL  CONSTANTS  (Y/N)  ?  *,\) 

BB  AD ( * , 30 )  RBSP 

C 

IF  (RESP  . BQ.  *  Y*  )  THEN 
C 

WRITE! *, 10) 

10  FORMAT!/, 2X,  ‘INPUT  THE  VALUBS  USING  DECIMAL  NOTATION 
DO  12  1=1,3 

WRITE (*,20)  LGBOM ( I ) 

READ!*, 70)  GBOM ( I ) 

12  CONTI  NUB 
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LBL(8)=’DBFL8CTION  AT  B  !  * 

C 

OPBN  ( 14 , FILB= ’ BCLBL. DAT'  , STATUS=  * NBW  ‘ } 

DO  155  1=1,8 

WRITB( 14, 156)  LBL(I) 

155  CONTINUE 

156  FORMAT ( A22  > 

CLOSE (  14) 

C 

WRITE ( * , 160) 

160  FORMAT( //, 6X , ’WOULD  YOU  LIKE  TO  SBE  THE  CURRENT  PLATE  BOUNDARY 
-  ,/,2X, 'CONDITIONS  ( Y/N )  ?  ' \ ) 

RE AD ( * , 30 )  RBSP 
C 

IF  ( RBSP  .80.  * Y* )  THBN 

OPBN ( 11 , FILB=’BC. DAT’ , STATUS = ’OLD’ ) 

OPBN ( 12, FI  LB = ’BCPOS.DAT’  , STATUS =’ OLD ’  ) 

WR IT8( * , 30)  *  * 

DO  165  1=1,4 

RBAD(11,70)  BC(I) 

165  CONTINUE 

DO  190  1=1,8 

R£AD{ 12,100  BCPOS(I) 

IF  ( BC POS (  I)  .GT.  0)  THEN 

WR I TB ( * , 186)  LBL(I) ,BC(BCPOS( I ) ) 

186  FORMAT(2X, A22, F20.  10) 

RND1F 

190  CONTINUE 

C  LOSE ( 1 1 i 
CLOSE ( 12 
END  I  * 

C 

WR I TE ( * , 17) 

17  FORMAT '//, 6X , ’ DO  YOU  WISH  TO  INPUT  OR  CHANGE  PLATE  BOUNDARY’, 

♦  /, 2X, ’CONDITIONS  (Y/N)  ?  ’,\) 

READ!*, 30:  RBSP 
C 

IF  ( RESP  . EO.  *  Y ’  )  THEN 
WR  I  TB ( » ,  19) 

19  FORMAT( //, 6X INDICATE  WHETHER  THE  BOUNDARY  CONDITIONS  FOR  THE 
$  / , 2 X  ,  ’PLATE  ARE  KNOWN  OR  UNKNOWN  BY  TYPING  IN  THB  LBTTBR  "U" 
S  / , 2 X ,  ’  FOR  UNKNOWN  BOUNDARY  CONDITIONS  AND  HITTING  RETURN  FOR 
W  /.2X, ’KNOWN  BOUNDARY  CONDITIONS’,//) 

C 

WR I TB ( *,20)  ’ SHEAR  AT  A  ? 

READ ( * , 30 )  8FIAGH  I  j 
WRITB ( * , 20 )  ’MOMENT  AT  A  ’ 

RE  AD ( * , 30  )  B  FLAG  1(21 

WR I TE ( * , 20  )  'SLOPE  AT  A  ” 

RE  AD  (  *  ,  30  )  BFLACK3) 

WR I TB !  * , 2  0)  ’ DBFLBCTTON  AT  A  ^ 

RE  AD  (  *  ,  30 )  B  F  LAG  1 ( 4  ) 

C 

WRITE (*,20)  ’SHBAR  AT  B  ^ 

READ!  * , 30  i  BFLAG2!  1 ; 
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r 


WRITB(*,20)  * MOMBNT  AT  B  ’ 

RBAD!*,30)  8  FLAG2 ( 2 ) 

WB I TB  f  * , 20  5  ’SLOPE  AT  B  ’ 

RB AD ( * , 30 )  8 FL AG2 ( 3  > 

WRITE!*, 20)  ’DEFLECTION  AT  B  ?  ’ 

READ!*. 30)  BFLAG2M) 

20  FORMAT ( 2X , A22 , \ i 
30  FORMAT !  A1  ) 

C 

WRITE ! * , 40 ) 

40  FORMAT://, 6X, ’TYPE  IN  THE  KNOWN  BOUNDARY  CONDITIONS  USING 
@2X , ' DECIMAL  NOTATION  WHEN  PROMPTED ’ , // ) 

C 

OPEN \ 8 , FILB= ' BFLaGI .DAT’ , STATUS = ’NEW ) 

OPBN { 9 , FI  LB*  * BFLAG2 . DAT’ , S T ATUS = ’ NEW '  ) 

OPEN <10 , FI LE= ’ BC . DAT'  . STATUS = 'NEW'  i 

C 

DO  50  1-1,4 

WR ITE ( 8 , 30 }  BFLAGI! I) 

IF  ! E  FLAG  HI)  .NE.  ’U’)  THEN 
L  =  L+  1 

WB I TE ( * , 20 )  LBLCI) 

RB AD ( * , 70 )  BC { L ) 

WB I TB (10,70)  BC ( L  > 

BNDIF 

50  CONTINUE 
C 

DO  60  1=1,4 

WRITS <  9, 30 )  BFLAC2 IJ) 

IF  ( BFLAG2 ( I )  . NE.  ’ U  * )  THEN 
L=  L*  I 

WRITE!*, 20)  LB  L  1*4) 

READ!*, 70)  BC(L) 

WRITE! 10,70)  BC ( L ) 

ENDIF 

60  CONTINUE 
C 

70  FORMAT! F20. 10) 

C 

CLOSE ( 8  ) 

CLOSE ! 9 ) 

CLOSE! 10) 

C 

OPBN! 1 1 , F I  LB  = ’ POS . DAT’  , STATUS = ’ NEW’  ) 

OPEN! 12, F I LE  = ’ BCPOS . DAT’  , STATUS = 'NEW* } 

C 

L--0 

DO  80  1=1.4 

IF  (BFLAGI ( I •  .EQ.  ’I”  THEN 
L  =  L  *  1 
POS i  I  •  =  L 
BCPos-,:  =0 
ELSE 
M  =  M-  1 

BCPOS  T ' = M 


98 


c 


c 

80 

c 


c 


c 

90 


100 

C 

C* 

c  *  *  * 
c* 

Cl 

ci 

c* 

C* 

Cl 

Cl 

c* 

Cl 

Cl 

Cl 

Cl 

Cl 

Cl 

Cl 

Cl 

c* 

Cl 

c* 

Cl 

Cl 

CI 

Cl 

c* 


POS( I  }=0 
BNDIF 

write (ii, 100'  pos( r : 

WRITE i 12 . 100 '  BCPOS / I . 

CONTINUE 
DO  90  1-5,8 

IF  ( BFLAG2 l I -4 i  -BO.  'U*>  THEN 
L=  L-  1 
POS (I!SL 
BCPOS ( I  ■  =  C 
ELSE 
M  =  M*  1 

BCPOS 1  I ) =M 
POS(  I  ) =0 
BNDIF 

WRITE( 11 , 100 )  POS(I) 

WRITE( 12. 100)  BCPOS { I ) 

CONTINUE 
CLOSE ( 11) 

C  LOSE (12) 

FORMAT*  ID 
BNDIF 


t*l**S«*S**ttll!lllll*****ll***ll*l**l*ll*!*ltll*ll*llltllltll» 

INPUT  CALCULATION  ADJUSTMENT  ITEMS.  N  IS  THE  NUMBER  OF 
INPUT  PARAMETERS  IN  THIS  SUBMODULE.  OBSERVE  THAT  SUBSEQUENT  TO 
THE  DB VBLOPMENT  OF  THIS  MODULE,  PARAMETER  NUMBER  3  WAS  NO  LONGER 
ACCESSED  FROM  ANY  OTHER  PROGRAM,  AND  THEREFORE  SERVES  ONLY  AS  A 
DUMMY  PARAMETER.  THE  REMAINING  PARAMETERS  HAVE  THE  FOLLOWING 
MBANINGS 

1.  EPSILON  REFERS  TO  THE  DISTANCE  BETWBEN  THE  RING  LOAD 
AND  THE  BOUNDARY.  NOTE  THAT  IN  SOME  INSTANCES,  LBTTING 
EPSILON  EQUAL  TO  ZERO  CAUSES  A  DIVIDE  BY  ZERO  TYPE  ERROR. 

2.  THE  INCREMENTS  FOR  R  REFER  TO  THE  STATIONS  ACROSS 
THB  PLATB.  THIS  APPLIES  TO  THE  NUMBER  OF  INTEGRATION 
INCREMENTS  FOR  DISTRIBUTED  LOADS,  AS  WELL  AS  THE  NUMBER 
OF  DATA  POINTS  FOR  THE  RESULTS. 

3.  THIS  PARAMETER  IS  NO  LONGER  USED 

4.  THIS  REFERS  TO  THE  EXTENT  BY  WHICH  THE  GRFENS 
riSCTlOS  PLATE  GEOMETRY  OVERHANGS  THE  ACTUAL  PLATE. 

FOR  EXAMPLE.  A  VALUE  OF  0.5  MEANS  THAT  THE  INNER  RADIUS  OF 
THF  RISC  LOADED  PLATE  DESCRIBED  BY  THE  GREEN'S  FUNCTION 
OVERHANGS  THE  ACTUAL  PLATE  BY  0.5  UNITS 
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c* 
c* 
c* 
c» 
c* 
c» 
c* 
c« 
c* 

Cl 
c* 

c» 

c 

N  =  ? 

LAD  J  ( 1 )  = * EPSILON  FOR  RO  =  ’ 

LADJ(2)= * INCREMENTS  FOR  R  =  * 

LADJ ( 3) - ’WIDTH  INTG  INCR  =  * 

LADJ(4)=*B0UNDS  FOR  A  AND  B  *  * 

LADJ(5)= ’ load  INCREMENT  =  ' 

LADJ(6)=’DBPTH  INTG  INCH  =  1 
LADJ(7)=’X  RMS  CONVERGENCE  =  * 

C 

WRITE}*. 102) 

102  FORMAT} //, 6X, ’WOULD  YOU  LIKE  TO  SEE  THE  CURRENT  ADJUSTABLE' 

*  ,/,  2X ,  ‘  parameter  settings'*  » • 

READ  t * , 30  >  RESP 

IF  < RESP  .80.  *  Y • )  THEN 

OPEN' 16,'FILE  = ’ ADJ . DAT’  , STATUS  * ’ OLD ’  ) 

WRITE}*, 30  '  ' 

DO  108  1=1, N 

READ ( 16,70)  ADJ(I) 

WRITE}*, 106)  LADJ' I),ADJ(I) 

106  FORMATi2X, A25, F20. 10 ; 

108  CONTINUE 

CLOSE} 16) 

END  I  F 
C 

WRITE} 1,110) 

110  FORMAT} // ,6X. 'DO  YOU  WISH  TO  INPUT  OR  CHANGE  ANY  COMPUTAT I ON A L ’ , 

♦  /, 2X,  1  ADJUSTMENT  ITEMS  (Y/Ni  ?  ’ ,\) 

READd, 30  RESP 

C 

IF  RESP  . EO.  ’ Y*  THEN 
WRITE  *,120 

120  FORMAT  INPUT  THE  VALUES  USING  DECIMAL  NOTATION 

DO  125  I  -  1  .  S 

WRITE  * .  126  IADJ .  T  • 

READ  *  .70  A DJ  I  . 

: 25  CONTINUE 


5.  THE  LOAD  INCREMENT  IS  THE  PERCENTAGE  OF  THE 
ORIGINALLY  SPECIFIED  LOAD  MAGNITUDE  THAT  IS  APPLIBD  DURING 
EACH  LOAD  INCREMENT  AS  THE  PLATE  IS  LOADED  IN  THE  PLASTIC 
RANGE 

6.  raiS  IS  TBE  NUMBER  OF  INTEGRATION  INCREMENTS  USED  TO 
DETERMINE  TBE  PLASTIC  MOMENT  ACROSS  THE  PLATE. 

7.  RMS  IS  THE  ROOT  MEAN  SQUARE  OF  THE  PLASTIC  MOMENT 
INCREMENTS  FOR  EVERY  POINT  IN  THE  PLATE  WHERE  THE  PLATE  IS 
NO  LONGER  ELASTIC.  THIS  IS  COMPARED  WITH  THE  RMS  VALUE  FOR 
THE  PREVIOUS  INCREMENT,  TO  SEE  WHETHER  THE  SOLUTION  HAS 
SUFFICIENTLY  CONVERGED.  1*  CONVERGENCE  MEANS  THESE  NUMBERS 
AGREE  TO  WITHIN  LESS  THAN  IX 
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126 


FORMAT; 2X , A25 


OPEN ; 13 , FILE*  * ADJ . DAT* , STATUS- ’NEW’ ) 

DO  130  1*1. N 

WRITE! 13,70)  ADJ(I) 

CONTINUE 
CLOSE {13} 

END  I  F 

WR I TB ( *  135/ 

format; ///, 1  YOU  ARE  LEAVING  TRE  INPUT  MODULE'.///) 
END 
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c 

<:****>****************************** ts***«t*t*t*s*** ******************* 

c* 

C*  PROGRAM  LOAD 

C* 

C*  THIS  PROGRAM  ALLOWS  THE  USER  TO  INPUT  THE  DESIRED  LOADING 

C*  CONDITION  FOR  THE  PLATE.  UNIFORM  DISTRIBUTED  LOADS,  RAMPS,  AND 
C*  PARABOLICALLY  DISTRIBUTED  LOADS  ARE  PERMITTED.  AS  WELL  AS  RING 
C*  LOADS.  TO  RUN,  THIS  PROGRAM  MUST  ACCESS  THE  FILE  LOAD.DAT. 

C* 

C*  USING  THE  INPUT  DATA,  THE  PROGRAM  CALCULATES  THE  TOTAL  LOAD 

C*  PER  UNIT  LENGTH  THAT  ACTS  ON  A  RING  OF  THE  PLATE  SURFACE  OP  A 
C*  WIDTH  DETERMINED  BY  ADJUSTABLE  PARAMETER  #3  OF  THE  INPUT 
C*  PROGRAM 

C* 

C*s*************t*t**t***********i************t*t********************* 

C* 

c 

PROGRAM  LOAD 
C 

IMPLICIT  REAL*8  (A-H.0-2) 

CHAR ACTER* 1  R8SP 
CHARACTER*30  L B L { 4 ) 

REAL*8  LD  f  50 } 

INTEGER  N 
C 

LBL; I )  = ’MAGNITUDE  AT  OUTER  RADIUS  ==  ’ 

LBL(2} =’ DISTANCE  FROM  PLATE  CENTER  =  ' 

LBL  3 )= 'MAGNITUDE  AT  INNER  RADIUS  =  * 

LBL(4) =’ DISTANCE  FROM  PLATE  CENTER  =  * 

C 

OPEN, 6, FI  LB  = *  LOAD . DAT* , STATUS = ’OLD’  } 

RE  AD  (  6 , 90  )  LD 
N=INT(6.0-2.0*LD(6) ) 

CLOSE (6' 

C 

WRITE <  * , 10  j 

10  FORMAT ( /  / ,  ’  DO  YOU  WISH  TO  SEE  THE  LOADING  CONDITIONS'. 

-  /,’  (Y/N)  ?  \\) 

READ (*,60;  RESP 
IF  ( RBSP  . EQ.  ' Y’ )  THEN 
write; *, 15) 

15  FORMAT!//,  ’  THE  EXISTING  LOADING  CONDITIONS  FOLLOW  i 

DO  40  1-1,4 

WR I TB ( * , 3  0 )  LB  L ( I ) , LD ( I J 
30  FORMAT ( 2X , A30 , F20  .  10  ' 

40  CONTINUE 

WRITE  *,6C  ’  ’ 

IF  . LD . 5  . EO.  1.0'  THEN 

WRITE*  * , 80  'LOADING  IS  LINEAR' 

ELSE 

WRITE  *,80  'LOADING  IS  PARABOLIC ’ 

END  I  F 

WRITE  *.60 

J  -  5 
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43 

45 


C 


46 


C 


50 


60 

C 


DO  45  J  =  ? , 6*  INT { LD ( 6) ) 

J  =  J  +  2 

WH1T£[®*60'  *  * 

WRITE  <  *  [  43  '  ’MAGNITUDE  Of  RING  LOAD‘,1-6, 
WRITS(*,43)  ’LOCATION  OF  RING  LOAD ’  ,  I  -  6 , 
FORMAT (2X,A24,I3,A3,F20. 10: 

CONTINUE 
END  I F 


’  ,  LD  ■  J 
, L V:J~ 


FORMAT* /!6:x. ’DO  YOU  WISH  TO  CHANGE  THE  LOADING  CONDITIONS 
*  /, * '  ( Y/N  )  ?  ’ , \) 

rsao<»,6o;  resp 


IF  (BESP  .BO.  ’Y’l  THEN 

F08M*T(//?2*. ’00  YOU  WISH  10  CHANGE  THE  DISTRIBUTED  LOAD'. 
♦  /,’  (Y/N)  •>  ’\) 

BEAD(*,60>  HESP 
FOHMAT(Al) 

IF  (HESP  BO.  ’Y’)  THEN 


70  FOMAt!//°2X.  'SPECIFY  THE  LOAD  PROFILE  USING  DECIMAL', 

.  /  _  •  NOTATION  ... 

DO  75  1=1,4 

tfRITEi».80)  LSUl; 

BEAD' » , 90  LD' I  1 
75  CONTINUE 

80  FORMAT i 2X , A90 , \ • 

90  FORMATCFSO. 101 

C 

UD T (  *  J  00  ) 

100  FORMAT;// ,2*. ■ IS  THE  load  linear  or  PARABOLIC  ll/P 
READ ( * , 60  )  RESP 
IF  (RESP  . EQ .  ' L ’ )  THEN 
LDi51  =  1-0 
ELSE 

10(5)  -  2.0 

END  I  F 
END  I F 

C 

110  FORMAT*  z/!?* . ' DO  YOU  WISH  TO  MODIFY  ANY  RING  LOADS  (Y  N 
READ  * , 60  )  RESP 
I f  i RESP  .  EQ.  * Y’  THEN 

wfi’TF ■ •  120  , 

120  FORMAT: '  ,2X,’THE  NUMBER  OF  RING  LOaDS  (INTEGER 

READ  * , 130  N 
130  FORMAT: 13 

LD  6  -  REAL : N 
K  5 

00  20 0  1  r  1  ,N 
K  K  *  2 
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WRITE  » , 60  ’  * 

WRITE  *,210  'MAGNITUDE  OF  BING  LOAD '  .  I  .  ’a 
READ  •. 90  LD  K 

WRITE  *,210  ’LOCATION  OF  RING  LOAD  ’.I.  ’  = 
READ  * , 90 :  LD  K*1 } 

0  CONTINUE 

C  FORMAT  2 X . A23 . 13 , A2 , '  ) 

END  1  F 


OPEN , 7, FILE:  ’LOAD. DAT’  , STATUS  r ’ NEW  ) 

DO  219  I M . 50 
WRITE  7,90  LD  I 
9  CONTINUE 
END  I  F 

WRITE  *,220 

0  FORMAT,  6X ,' YOU  ARE  LEAVING  THE  LOADING  MODULE’ 
END 
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c 

c* 

CfllltttlttflitttMUttttltttUMtt 


c* 

C*  CALL  THE  SUBROUTINE  YLD .  WHICH  WILL  RETURN  A  LOAD  MATRIX 

C*  SUCH  THAT  YIBLDING  HAS  JUST  STARTED  TO  OCCUR  IN  THE  PLATE 


C* 

C«xx*axxxxxxxxxtxt*x*t*tt*«*t******sxxxxxxxxxxxx***x*xt*xttxxtxxttxi 


C« 


C 


CALL  YLD BE  , IMAX , LDMATI ; 


C 

C 


******  «mmmn 


REPEAT  THE  PROCESS,  AND  THEN  CALL  UPDATE  TO  REVISE  OUR  STATE 
MATRIX  TO  REFLECT  THE  CONDITION  OF  THE  PLATE  AT  THE  ONSET  OF 
YIBLDING 


C 


C 

C 


CALL  AMATR(BFLAG) , BFLAG2 , BC , BCPOS , POS , LDMATI , LHS , AMAT ; 
CALL  GAUSS' LHS, AMAT, X) 

CALL  BC»(POS, BC, X. BCMAT} 

CALL  SOLVER* LDMATI , BCMAT ) 

CALL  LPDATB(STT2 , IMAX. OUT, LDMATI , LL, LDMATO  ' 

CALL  RSSUKM.STTl.COUNT,  LDMATO 


NOW,  INPUT  AN  INCREMENTAL  LOAD,  AND  COMPUTE  THE  RESULTING 
MOMENTS , DEFLECTIONS  ETC,  USING  AS  AN  INPUT  THE  PLASTIC  MOMENTS 
(IF  ANY,  AND  OTHER  DATA  FROM  THE  PREVIOUS  LOAD  CASE 


C* 

C 

DO  no  1=1  ,  N N 

LDMATI  I)=ADJ(5  j  * LDMATO (  I  , 

110  CONTINUE 
C 

116  CONTINUE 
C 

C  WRITE*  *,  130 i  STT.  IMAX, 9 i ,STT; IMAX.  10 i 

C  130  FORMAT* 5X , ’ MRP  IN  r  ’ .F20.10,’  MTP  IN  =  ’.F20.10 


CALL  AMATR  B FI AG  1 . B FL AG2 . BC . BCPOS , POS , LDV AT  1 . LHS . AMA 7 
CALL  GAUSS  LHS , AMAT. X 
CALL  BCM  POS , BC , X . BCMAT 
CALL  SOLVER  LDMATI. BCMAT 
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c* 

c* 

c* 

c* 

c* 

c* 

c* 


GO  INTO  THE  SUBROUTINE  FLAST  TO  CALCULA  TE  THE  1NCREHESTM  PLASTK 
MOMENTS  AT  BACH  STATION  ACROSS  THE  PLATE  WHERE  YIELDING  HAS 
OCCURRED.  THE  PLASTIC  MOMENTS  THAT  ARE  RETURNED  BY  PLAST  ARE 
SQUARED  AND  ADDED  TO  THE  TOTAL  IN  THE  VARIABLE  RMS .  RMS  IS 
COMPARBD  TO  THE  RMS  VALUE  FROM  THE  PREVIOUS  ITERATION  (SAVED  AS 
RMS  i  i 


C * 

C 


200 

C 


RMS  1  =  RMS 
RMS  s  0 . 0 
IZ-NINT  .  AD J  •  6  ;• 

DO  200  I  =  1  , NN 

STT1 (  1 ,  1 i =  STT  ( 1,9. 

STT1 i l , 2 ) =STT( 1,10: 

CALL  PLAST ( I  ■ 

ir  (sbtli.izj  .eg.  geom(6M  tbbn 

RMS  =  RMS*STT 1  I , 9  j  *  * 2 . 0*STT  < 1,10  **2.0 
END  I  F 
CONTINUE 


WRITE  *,810  RMS 
810  FORMAT  5X RMS  =  ’.F20.10 
C 

C* 

C •  THIS  WAS  THROWS  IN  TO  AVOID  A  DIVIDE  BY  ZERO,  UNTIL  THE  PROCESS 

C*  HAS  BEEN  "PRIMED” 


C$ 

C 

IF  RMS  .80.  0.0:  THEN 
RMS  -  1  . 0 
E  NO  I  F 
C 

IF  (ABSSQRT(RMS)-SQRT(RMS1))/SQRT( RMS itlOO.O  .GT.  ADJC7I1  THEN 
DO  815  1=1, NN 

STTI I , 9 )  =  ( S TT ' 1 , 9 ) *STTl ( I,  1) )/2.0 
STT'  I  .  10  )*(  STT(  1 ,  10  WSTTK  1,2;  >  /  2 . 0 
815  CONTINUE 
GOTO  116 
ENDI  F 
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c* 

C«  CALL  THE  SUBROUTINE  UPDATE.  WHICH  WILL  UPDATE  OUB  STATE 

C*  MATRIX  AND  OUR  LOADING  MATRIX  TO  REFLECT  THE  CONDITION  OF  THE 
C*  PLATE  AT  THE  EN:  OF  THE  LOADING  INCREMENT.  THE  PARAMETER 
Ct  COUNT  KEEPS  TRA  JK  OF  LOAD  INCREMENTS  FOR  PURPOSES  OF  DISPLAYING 
C*  RESULTS.  WHEN  COUNT=10,  RESULTS  ARE  DISPLAYED  AND  WRITTEN  TO 
C*  FILE.  IE  BVERY  10  LOAD  INCREMENTS- 
C* 

C* 

C 

CALL  UPDATE  ; STT1 ,  IMAX  .  OUT,  1DMAT1  ,  J  L  ,  LDMATO 
CALL  RESUL.M. STT1 .COUNT. LDMATO 
IF  COUNT  EQ.  10  0  THEN 
COUNT =0 . 0 
END  IF 

COUNT  =  C0UNT- 1  0 
M- M*  1 

WRITE f  *,117.  M 

117  FORMAT; 5X, ’LOAD  INCR  *  ’,13? 

C 

C* 


************ 


C* 

C*  CHECK  TO  SEE  IF  EXCEEDING  ULTIMATE  STRAIN.  IF  SO,  KICK 

C*  OUT  OF  THE  PROCEDURE.  IF  NOT,  WE  SHOULD  INCREMENT  ON  UP  AGAIN 
BY  USING  AS  OUR  INPUT  THE  INCREMENTAL  LOAD  MATRIX. 


C 

C 

C 


************* 


IF  . OUT  . NE.  1 !  THEN 
GOTO  116 
END  IF 

CALL  RESULlM.STTl .COUNT. LDMATO 

CLOSE  15 

END 
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c 

c* 

c* 

C*  SUBROUTINE  INTEC 

C* 

C*  THIS  SUBROUTINE  IS  CALLED  BY  THE  PROGRAM  MAIN  TO  CREATE  THE 

C*  A8RAY  LDMAT.  THIS  ARRAY  CONTAINS  THE  VALUES  OF  THE  TOTAL  LOAD 
C*  TOR  EACH  SEGMENT  AlOSG  THE  PLATE  RADIUS.  IT  IS  CONSTRUCTED 

C»  BY  USING  THE  AVERAGE  OF  THE  TWO  ENDPOINTS  FOR  DISTRIBUTED 

C*  LOADS  AND  ADDING  TO  THIS  VALUE  THE  MAGNITUDE  OF  ANY  RING  LOADS 

C*  THAT  ARE  SITUATED  WITHIN  THE  SEGMENT. 

C* 


C* 

c 

SUBROUTINE  INTEC  LD, LDMAT 

C 

IMPLICIT  BEAL«8  :a-H,0-Z 
RE  A  L  * 8  LDMAT I  50 '  , ID i 50 } , NU 
INTEGER  N 

(INCLUDE : 'COMMON. FOR’ 

C 

A  =  GBOM( 1  - 
B  =  G80M ' 2 ' 

NU=GEOM< 3 
D = G ROM  4 
C 

K-N1NT  ADJ  2 
DEL:  A-B  aDJ  2 
R  -  B 

N 1 N  INT  ■  6 . 0-  2 . 0*  LD  '■  6  t  ' 

C 

C* 


C*  CHECK  TO  SEE  IF  THERE  IS  A  DISTRIBUTED  LOAD.  IF  NOT,  PROCEED 
C*  TO  THE  S8CTION  FOR  RING  LOADS 

C* 


C* 

C 

IF  LD1.  • NE .  0.0  .OR-  LD • 3 .  . NE .  0.0  THEN 
C 

R  -  B 

SLP-  LC  }  - LD  J  LD  2  • LD  4 
DO  IOC  t- 1 ,K 
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SECTION  FOR  LINEARLY  DISTRIBUTE!;  LOADS 


IF  LD  5  . EG.  I  ?  THEN 


8- R-DEL 


IF  R  LD  2  . GT.  DEL. 10 . 0  THFN 

LDMAT  I  =  0-L 

ELSE  !F  R  .GT.  LD  4  THEN 

LDHAT  I  -LP«  R-LD  4  -DEL. 2.0  'LD  3  ;*DE 

ELSE 

LDMATl I)=0.0 
ENDI  F 


itnlfittiiiinmtittMMtttu 

SECTION  FOR  PARABOLIC  LOADS 


ELSE  IF  ' LD- 5  . EO  2.0  i  THE* 

R  R-  DE  L 

P -  L D I  2  •  - L D ' 4  i  *#2.0  4.0/  LD  1  - LD ■  3 
!  F  .  LD  1  I  .  GT  LD-'  3  >  >  THEN 

Y= ' R-DEL/2. 0-LD; 4  **2.0/4. 0  P-LP  3 

ELSE 

Y= - • LD ■ 2 l-R-DEL/2 • 0 ) **2 . 0  4.0  P*LD  4 
END  IF 

IF  R  -  L  D  (  2  '  -GT.  DBL'10.0)  THEN 
LDMAT'  I  ) =0 . 0 

ELSE  ! F  ■ R  . GT .  LD  4  THES 
LDMAT  I  -  Y  *  D F.  L 

ELSE 

LDMAT  1  -0.0 

END  I  F 
EM  l  F 


■  1 1  u  1 1 1  > 


1  1  1 


r 

c 

1 00  CONTINUE 
END  I  F 


J  =  5 

IF  N  .NE.  0)  THEN 


|  DO  300  I  =  7 , S-N INTi LD • 6  * 

;  J  =  Jf  2 

I  L  =  I N  T  i  ( LD  ?  J  *  1  »  -  P  '»  /  D  B  L  J  ♦  1 

i  LDMATl  Li  =  LDMAT(  L)-*-LD(  J  ) 

'  300  CONTINUE 

END  I  F 

I  i 

RETURN 

END 


C* 

C* 

Zt 

C* 

C* 

C* 

C* 

C*t 

ct 

c 


SVBJfOU?  1NB  AMATR 

THIS  SUBROUTINE  IS  CALLED  BY  THE  PROGRAM  MAIN  TO  CREATE  THE 
4X4  A  MATRIX.  WHICH  IS  SOLVED  USING  GAUSS  ELIMINATION  TO 
DETERMINE  THE  UNKNOWN  BOUNDARY  CONDITIONS. 


SUBROUTINE  AMATR ; B  FLAG  1 , BFLAG2 , BC. BCPOS , POS , LDMAT  3  ,  IHS . AM AT ! 


IMPLICIT  REAL*8  A-H, 
CHARACTER* I  BFLACl  4 
R E  A L * 8  BC . 4  ,  LHS  4  .  : 
AMAT  4,4  .LDMAT 
*  STTAV  50.2 

INTEGER  POS  8 ■ . BCPOS : 

S INCLUDE. ’COMMON. FOR’ 

C 


O-Z 

,  BFLAG2  4 

NC.NC.GRNMT 1  4.6  , GRNM72 . 4  .  e  . 
4  , LDTCT  4  , LDMATl'50  , M F R T 0 T 

8  .  M  .  N 


A-GEOM  3 
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B  =  G  BOM ( 2 ; 

NU-GBOM i 3 ' 
d-gbom;4 , 

INC=ADJ( 1  ' 

c 

DO  101  1=1.4 
LHS l I ) =0 . 0 

rr  (I  .  gQ.  1  .OR.  I  .80.  3 1  THEN 
RO= A- INC 
ELSE 

BO=B+ INC 
END  I F 

C 

H  =  A 

CALL  GRNFNC l R, RO, NU, D , GRNMTi 
R  =  B 

CALL  GRNFNC ( 8. RO , NU , D , GRNMT2 
C 

ct 


C«  REDUCE  THE  TBRMS  ORIGINALLY  ON  THE  RBS  TO  A  4X4  MATRIX 

C*  COEFFICIENTS  (AMAT;  AND  A  4X1  MATRIX  OF  NON-COBFf 1 C I  ENT 
C*  TERMS  (  LHS  ) 

C* 

C* 

C 

IF  (I  .  EQ.  1  .OR.  I  .80.  2.)  THEN 
K  =  0 
L=  0 

DO  90  J  =  1 , 4 

IF  (BFLAGliJ)  .  NE.  ’I")  THEN 
K  =  K-  1 

LHS  i  I  =GBNMT1  :  J  ,  1  ;  *BC  K  '  *  A  *  •  -  1  •  0  **  J-li-LHSil) 

ELSE 

L=L-1 

AMAT ( I , L ) =GBNMT 1  (  J  , 1)»A«(-I.0!MU) 

8NDIF 

73  F0RMAT'2X, F10. 5) 

90  CONTINUE 

C 

DO  100  J= 1 , 4 

IF  t BFLAG2 ' J )  . NK .  ’U’  '  THEN 

E-  S*  I 

LHS  ]  •  sGRNMT2 .  J  .  1  *BC  K  »B»  -)  •  *J*!.HS  I 
ELSE 

L  -  L  -  1 

AMAT  I  , L  =  G  H  N  M  T  2  >  J  .  1  *  B  *  1  •*  J-i 

BSD  IF 

100  CONTINUE 


ELSr 
K  :  0 
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L  =  0 

DO  no  J=  1 , 4 

IF  (BFLAGliJ)  . NE .  *0’)  THEN 

E  =  S*l 

LHSv'I.'=Gfi.VMTl  { J,2)*BC(K)*A*(-1 . 0/**{J-l  }-LR$(  I 
ELSE 

L=  L-*- 1 

AMAT  ( I , L ) -GRNHT l(Jf2)*A*(-1.0)**(J) 

END  I  F 

110  CONTINUE 

C 

DO  115  J= 1 , 4 

IF  ( B  FLAG2 ( J  )  -NE.  * U 1 )  THEN 
1  =  1+1 

LBS ( I ) sGRNMT2 i J , 2 } *BC { l ) *B * ( - 1 ) ** ( J } ♦ LHS ( I } 
ELSE 
L=  L*  1 

AM AT ( I , L ) =  GRNMT2 ( J , 2 • *B  * ( - 1 ) #  * ( J -  1  ) 

END  I  F 

115  CONTINUE 

END  I  F 


C 

C* 


MOVE  THE  TBRMS  ORIGINALLY  ON  THE  LBS  TO  THE  RHS ,  FORM  A  NEW  AMAT 
MATRIX  AND  A  NEW  LHS  MATRIX  WHICH  CORRESPOND  TO  THE  A  AND  B  MATRICE 
OF  THE  SYSTEM  AX=B 


C* 

IF  . B  F  L  A  G 1 \ 4  ;  .NE.  ’I”  .AND.  I  . EQ .  1)  THEN 

LHS  (  1  >=LBSi  1 )  «- RO •  B C  v  BCPOS. 4  ;  ; 

ELSE  IF  >  BFLAG1 ( 4 )  .EQ.  ’U’  .AND.  I  . EQ .  1)  THEN 

AMAT (  I , POS ( 4 >  ) =  AMAT (  I,P0S(4'  )-R0 

END  IF 
C 

IF  ( BFLAG2 ( 4 )  .NE.  ’IT  -AND.  I  . EQ .  2)  THEN 
LHS ( *  ■  =  LHS { 2 ) *RO*BC ( BCPOS ( 8 ) ; 

ELSE  IF  ( B  F  LAG2 ( 4 )  .EQ.  ’U*  .AND.  I  .EQ.  2)  THEN 
AMAT( I, POS (8) )=AMATt. I , POS ( 8 ) J-RO 

BNDIF 

C 

IF  (BFLAGK3)  .NE.  ’IT  .AND.  I  .  EQ .  3>  THEN 
LHS ( 3 ) - LHS ( 3 ) - RO*BC < BCPOS ( 3 ) } 

END  IF 

IF  BFLAG1  4!  . NF .  ’U’  .AND.  I  . EQ .  3  THEN 

LHS  3  =  LHS  3  -BC  BCPOS  4 

END  I  f 

IF  B  FLAG;.  3.  .  EQ .  T*  .  AND.  I  .EQ.  3  THEN 

AMAT  I , POS .  3 '  =  AMAT .  I . POS . 3 ;  ‘RO 

END  I  F 

IF  B  F  LAG  1  4  .EQ.  '  V  *  .AND.  I  .EQ.  3  THEN 

AMAT  I  . POS ■  4  =  AMAT  I  . POS  4  1.0 

END  I  F 
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f 


IF  ( BFLAG2 ( 3 )  .NB.  ’ U ’  .AND.  I 
LHS ; 4  >  =  IBS ' 4 } -RO«BC { BCPOS { 7 ) ) 
BNDIF 

.  BO. 

4  1 

THEN 

IF  ( B FLAG2 ( 4 )  . NB .  'U'  .AND.  I 

L8S (4>- LHS ( 4 ) +BC ( BCPOS  (  8  1  - 
END  I  F 

.  BO. 

4) 

THBN 

IF  ( BFLAG2 ( 3 )  .BQ.  ’U'  .AND.  I 
AMAT( I , POS (7) )  =  AMAT ( I , POS  (  7 ) ) 
BNDIF 

.  BO. 
♦  RO 

4 

THEN 

IF  ( BFLAG2 * 4 ;  .BO.  *U’  .AND.  I 

.  BO. 

4 

TEEN 

AMAT( I , POS '8) )=AMAT( I , POS : 8: )-l . 0 
END  I  F 
C 

101  CONTINUE 
C 


C« 

C ***«*« **«*«>* *«$»******>«* *«t**t**t«t* *9******* *«**«*» ******** 

c* 

C*  FINALLY,  ADJUST  THE  LBS  TO  ACCOUNT  FOR  THE  EXTERNAL  LOAD 

C*  AND  PLASTICITY  BY  CALCULATING  AND  SUBTRACTING  APPROPRIATE 
C*  TBRMS . 

C* 


Kr  NINT ( AD J ( 2 )  ) 

DELMA-B  ,  /ADJ(2) 

C 

DO  500  1  =  1  .  K 

R  =  B-  DEL/2. 0 ♦ R  E  A  L ( I )  *  D  E  L 
HO  =  A- ADJ (  1  ) 

5TTAV ( I , 1 ) = ( STT ( I , 9 ) ♦STT{ I ♦ 1 ,911/2.0 
STTAV( 1,2)  =  f STT{ 1 , 101-STT( 1*1, 1011/2.0 
Call  GHNFNC(R,RO,NU, D.GRNMTl) 

LDTOT(  I )* LDTOT ( 1 ) * IDMAT1 ( I ) tGRVMTl ( 1 ,  1  } *R 
MPRTOT { 1 )=HPRT0T( 1)*STTAV{ 1 , 1  1  * ( -GRNMT 1 ( 3 , 1 ) / D ♦ 
NU/ R*GRMMT 1(2, 1) ) *0BL*R 

MPTTOT ( 1 }=MPTT0T( 1 j+STTAV( I ,2) *( -GRNMTl (2, 1 ) %DE  L  > 
LDTOT { 3 ) = LDTOT ( 3 ) ♦ LDMAT 1 ( I ) *GRNMT 1 ( 1 , 2 ; *R 
MPRTOT(3)=MPRTOT(3)*STTAV(  I  ,  1  )  *  {  -  G  R  NM  T 1  (  3,2  ,'  /  D  * 

♦  NU/R*GRNMT1  (2,2)  )*DEL*R 

MPTTOT* 3) =MPTTOT(3)-STTAV( I , 2 1  * ( -GRNMT 1 ( 2 , 2 ) *UEl ) 
RO  =  B  +  AD J (  1  ) 

CALL  GRNFNCf  R,  RO,  NU,  D,  GRNMTH 

LDTOT (  2  )  =  LDTOT  (  2  '  ♦LDMATl  (  I  )  *GRNMT  1  >'  1  ,  1' *R 


* 

MPRTOT i 2 ) -MPRTOT ' 2  *S  TTA  V  t 
NU / R*GRNMT 12,1'  *DE l»R 

I  .  3 

:•  •  :  -  GRNMT  1  '3,1 

D* 

MPTTOT  2  =MPTTOT. 2 1 -STTaV 
LDTOT  4  =  LDTOT , 4  -LDMAT1  1 

1.2*- GRNMT l  2,1 
•  GRNMT i  1,2  *  R 

•  DPI 

_ 

MPRTOT  4  =MPRTOT  4  *STTAV- 
NU/R*GRNMT1 ' 2,2 :  *DEL»R 

I  .  1 

»  -GRNMT '.3,2 

D  ■ 

MPTTOT. 4  '  =MPTTOT ( 4  *ST7AV: 

I  ,  2 

*  -GRNMT 1  2,2 

•  DEI 

500  CONTINUE 


C 


115 


► 


c 


DO  105  1=1,4 

LBS  (I  )  =  LHS(  I  )-LDTOT(  I  )  ♦NPRTOT  (  7  )  »-MP?TOT  '  I 
CONTINUE 

DO  1000  1=1,4 
LDTOT ( I ) =0 . 0 
HPRTOTC I)sO. 0 
MPTTOT' I) =0.0 
1000  CONTINUE 

RETURN 

END 


f 

i 


C 

c* 

Ct*»****t$t**tt*t*t*tt***t***«t*tttttt*sts***tis*9*ts*tt**t*«**tt* 

C* 

ct  SUBROUTINE  GRNFNC 

c* 

C»  THIS  SUBROUTINE  CALCULATES  THE  VALUES  OF  DEFLBCT I  ON ,  SLOPE, 

C*  MOMENT  AND  SHEAR  AND  ASSOCIATED  DERIVATIVES  FOR  THE  GREEN’S 
C*  FUNCTION,  AND  RBTURNS  THESE  VALUES  TO  TRE  MAIN  PROGRAM.  TO 
C*  INDICATE  DERIVATIVES  OF  PARAMETERS,  THE  NUMBER  1,  2.  OR  3  IS 
C*  TACKED  ON  TO  THE  END  OF  THE  VARIABLE  NAME.  HENCE.  THE  FIRST 
C*  DERIVATIVE  OF  L9  WRT  PO  IS  L91 
C* 

C» 

C***t**t***tt**%****tt*t«***t*t%ttt«t«ttt»**ttt*t»tt*t**t»«****tstst 

C* 

c 

SUBROUTINE  GRNFNC ( R , RO , NU , D . GRNMAT ) 

C 

IMPLICIT  REAL*8  (A-R.O'Z) 

RSAL*8  NU, GRNMAT (4, 6) , 13, L31, L32, L33, L9 , L9 l , L92 , L93 
$  JNCLL’DB  :  ’COMMON.  FOR’ 

C 

W=  l  .  0 

A=GBOM( 1 ) +  AD J ( 4  ) 

B  =  GBOMl 2 ) -  AD J ( 4 ) 

THB  =  0 . 0 
MB  0.0 
OB  -  0 . 0 
C 

IF  iR  .GT.  RO;  THEN 
RHO=  1  . 0 
ELSE 

RRQ  =  0 . 0 
end:  f 
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c 

c* 

c* 

C*  CALCULATE  PLATE  CONSTANTS  AND  PLATB  FUNCTIONS,  AND  THE  IB 

C*  DERIVATIVES 


C* 

C 


L3-RO/4.0/A*(  (  (RO/A)**2.  0+  1 . 0  )  *D  LOG  ■;  A  /  BO  )  +  (  RO  /  A  )  *  *  2 . 0- X  .  0  ) 

L31 «  1 . 0/4 . 0/ A*  (  D LOG  ( A/ RO } * l 3 . 0*RO* *2 . 0  /  A* *2 . 0*  1 . 0  )  *-2 . 0 * 

♦  f 80**2. 0/A**2 . 0-1 • 0 ) ) 

L32  =  1 . 0/4 . 0/ A* ( 6 . 0*RO/A**2 . 0*DLOG !  A  /  RO  }  ~  RO  /  A  *  *  2 .0-1 . 0 / BO ' 

L33s 1.0/4.  D/A*(60/A**2-0*: D LOG ( A / 80  -1 . 0 ■ *  1 . 0/A**2. 0*1.0 /RO 

*  **2.0) 

C 

L  9  =  RO / A  *  < ( 1 . 0 ♦ NU ) / 2 . 0*DLOG ( A/BO )  *  ( 1 ,0-NU)/4.0»(1.0-(RO/A)**2.0>  i 
L91 s l . 0/4 . 0/ A* <2. 09(1 . 0+NU ) *DLOG( A- BO/-(3.0*NU+1.0)-3.0*( 1 . 0-NU  < 

♦  *BO**2. 0/A**2. 0} 

L92  =  - { 1 . 0*NU) /2 . O/A/RO-3. 0*RO*( 1 . 0-NU > /2 . 0/A**3 . 0 
L93 - ( 1 . 0*NU)/2.0/A/R0**2. 0-3. 0*( 1 .0-NU)/2.0/A**3 . 0 
C 

ClM  1 .0*NU)/2.  0»B.'A*DLOG(  A/B  ,  M  1 . 0-NU1 /4 . 0*(  A/B-B/ A  ) 

C?s 1 . 0/2. 0*( 1 . 0-NU**2 . 0 ) *( A/B-B/A ) 

r 

G3-RO/4.0/R*(  (  (  RO  /  R  )  *  *2 . 0  *  1 . 0  )  *  D  LOG  (  R  /  RO  >  ■♦  (  RO  /  8  )  *  *2 . 0  -  1  .0)*RRO 
G31*l . 0/4. 0. R*f 3. 0*RO**2. 0/R**2. 0*DLOG:R/RO) *2. 0*RO* *2 . 0 / R* *2 . 0- 
-  D  LOG i R ' RO : -  2 . 0  1  *RRO 

C32M.04.0  R*  6.0tRO  R* *2  . 0*D  LOG  R  RO  *  RO  8**2 . 0-  1  .  0  -  RO  ■  *RRO 
G  3  3  =  1  0'4.0/R*:6.0.  R  *  *  2 . 0  *  .  D  LOG  R  RO  , -1  .  O':  *1 . 0- R*»2. 0*1 , 0'R0**2  .  0 

♦  *  RRO 
C 

G  6  -  RO  /  4 . 0  /  R  *  «'  {  RO  /  R  i  *  *  2 . 0  -  1 . 0  *  2 . 0  *  D  LOG  (  R  /  RO  /  )  *RRO 
G6l-1.0/4.0/R*i3.0»;RO**2.0/R**2.0-1.0 ; *2  - OtDLOGI R/RO  >  , *RRO 
G62  = ( 3 . 0*RO/2. 0/R**3. 0-1 . 0/2- 0 'HO/H ; *RRO 
G63=(3. 0/2. 0/R**3. 0*1.0/2.0/R/RO**2.0, *RRO 
C 

G9  =  RO/fi* ( (1.0*HV)/2. 0  *D  LOG ( R/HO ) ♦ ( 1 . 0  - NU ) / 4 . 0  * ( 1 . 0 - ( RO /R ) **2 . 0 ) ) 

♦  *RRO 

G91  =2.0/4.0/R*(2.0*< 1 . 0  +  NU) * i DLOG ( R, RO } - 1 . 0 ) ♦ (1 . 0-NU ) -3 . 0  * 

♦  ( 1 . 0-NU) *RO**2 . 0/R**2 . 0 > *RR0 

G92= -  1 . 0/4 . 0/R *<2.0$< 1 . 0+NU) /BO+6. 0*80  8**2. 0*(1. O-NU)) *HRO 
G93  =  l  .0/4. 0/R*C2. 0*(  1  0*NU)/R0*»2  0 -6 . 0 / R* * 2 . 0  * (  1 . 0 - NU ) ) * RRO 
C 

F) = ( 1 . 0*NU  /2.0*BR*DI.OG'R/B’*  1 . 0  -  HU )  '4.0*  R  B-B  R) 

F2  =  1  . 0  4 . 0* ,  1  . 0-  .  B  R  **2.0*.  1.0»2.0*ULOG-  R  B 

F3-B  4. OR*'  B  R  *  *  2  .  0  -  1 . 0  *  D  LOG  R  B  •*  B  R  **2.01.0 

F4-1  .0/2.0*;  i  10*N'l  i*B  R*  1.0-NL"«R  B 

F5= 1 . 0/2. 0*1.0*  B  R  * *2 . 0 

F6-B'4.0/R*'  B  R  **2 . 0- 1  . 0*2.0*DLOG  R  B 

F7- 1 . 0/2. 0* ' i • 0  NU**2 . 0  *  R  B-B  R 

F8M  0-2.0*  1  ,  O-Nl  *  l  1  .  0-NU  ■  *  B  R  **2  0 

F9-e  Rr  i.0*.vr  2  .  O«f/LOG  R  B  *.i.C-VL  -J  .  0*  1.0  B  R  **2.0 
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c 

c* 

c* 

C*  CALCULATE  TBE  DEFLECTION,  SLOPE,  MOMENT  AND  SHEAR,  AND 

C*  ASSOCIATED  DERIVATIVES,  AND  STORE  IN  THE  MATRIX  GRNMAT. 

C*  THIS  MATRIX  IS  A  4  BY  4  MATRIX,  WHERE  ROW  POSITIONS  ARB  AS 
C*  FOLLOWS 

C* 

C*  l  DEFLBCT1 ON 

C»  2  SLOPE 

C*  3  MOMENT 

C*  4  SHEAR 

C* 

C*  AND  COLUMN  POSITIONS  ARE  AS  FOLLOWS 

C* 

C*  1  GREEN’S  FUNCTION  (NO  DERIVATIVES) 

C*  2  1ST  DERIVATIVE  OF  GREEN’S  FUNCTION  WRT  RO 

C#  3  2ND  DERIVATIVE  OF  GREEN’S  FUNCTION  WRT  RO 

C*  4  3RD  DERIVATIVE  OF  GREEN’S  FUNCTION  WRT  RO 

C* 

Ct  FOR  EXAMPLE,  THB  VALUE  RETURNED  BY  GRNMAT (2,3)  REPRESENTS 

C*  THE  SECOND  DERIVATIVE  OF  THE  SLOPE  FOR  THE  GREEN’S  FUNCTION 
C*  WRT  RO 

C* 

C* 

C 

YB-  -w*A*«3. 0  '  D*  i  C  I  *  L9/C7  -  L3  '< 

THB=W*a**2 . 0/D  C7*L9 

GRNMAT < 1, I '=~'YB-THB*R*F]-W*R**3.0'D*G3 
GRNMAT ( 2 , 1 ) =THB»F4-W»R**2. 0, D*G6 
GRNMAT ( 3 , 1 > - THB * D / R« F7 - W *R*G9 
GRNMAT(4, 1 ) -W9R0,  R*RRO 
C 
C 

YB1=-W»A**3.0/D*(C1*L91/C7-L3n 
THB1 =W*A«*2 . 0/D/C7*L91 

GRNMAT ( 1 ,2) =-( YB 1 ♦ THB 1 *R* F 1  - W*RA *3 . 0 / D *G3 1 ) 

GRNMAT ( 2, 2) = THB 1  * F4 - W * R **2 . 0/D*G61 
GRNMAT { 3, 2) =THB 1*0  R*F7-W*R*G91 
GRNMAT(4,2)*W/R*RRO 
C 
C 

YB2=-W*A**3.0/D»<C1 /C7*L92- L32 
THB2-W»A**2. O/D  C7*L92 

GRNMAT  ••  1  ,  3  •  s  -  Y  B  2  -  R  *  F  1  *  THB  2  W  *  R  *  •  3 . 0  .  D  *  G  J  2 
GRNMAT l 2 , 3  = F4 * T HB 2 - W # R t *  2 . 0  D»G62 
GRNMAT .3,3  ‘  1 D»  F7  fi«THB2W#R*G92 
GRNMAT  4,3  :  0 
C 
C 

YB3 -  -W*A **3 . 0  D*  Ci  C7  »  L93 - L33 
THB3--Wt  A*t2  .  0  D-C?»L93 

GRNMAT:  1,4  --  YB3-R4F 1 *THH3 - Wifi* *2 . 0  D*GS3 
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CRNMAT(2 , 4)=F4*THB3-W*R**2. 0/D*G63 
GRNMAT(3, 4 )=D*F7/R«THB3-W*R*G93 
G8NMAT ( 4 , 4 ) =0 


RETURN 

END 


C 

C* 

c* 

Ct  SUBROUTINE  GAUSS 

Ct 

Ct  THIS  SUBROUTINE  USES  GAUSS  ELIMINATION  TO  DETERMINE  THE 

C*  UNKNOWN  BOUNDARY  CONDITIONS.  THE  MATRICES  A  AND  B  OF  THE 

Ct  EQUATION  A X  - B  ARB  PROVIDBD  BY  THE  SUBROUTINE  AMATR.  THIS 

C*  SUBROUTINE  THEN  RETURNS  THE  MATRIX  X  OF  UNKNOWN  BOUNDARY 
Ct  CONDITIONS.  THIS  SUBROUTINE  WAS  TAKEN  FROM  THE  BOOK 

C*  "NUMERICAL  ANALYSIS"  BY  L.W.  JOHNSON  AND  R.D.  RIESS. 

Ct 

Ct 

C 

SUBROUTINE  GAUSS, B, A, X) 

C 

IMPLICIT  REAL*8  ;a-H,0-Z 
CHARACTER*22  BCLBL-8) 

f?EAL*8  a;  4, 4  )  ,B  f  4.)  ,  X(  4  ,  BCMAT(  8  )  ,  BCi  4  ■ 

DIMENSION  AUG  v  50 , 51 
INTEGER  PO  S  (  8  ,  B  C  PO  S  8  '  ,  M  ,  N 
C 

NM 1  =  3 
NP  =  5 
N  -  4 
C 

C  SET  UP  THB  AUGMENTED  MATRIX  FOR  AX  =  B 
C 

00  2  1=2, N 
DO  1  J  =  1  , N 
AUG' I , J  >=A\ I . J ) 

3  CONTINUE 

AUG  I , NP  : B ■ I 
2  CONTINUE 
C 

P  THF  01  TER  LOOP  USES  ELEMENTARY  ROW  OPERATIONS  TO  TRANSFORM 

:  THE  augmented  matrix  to  echelon  form 
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c 

DO  8  1=] . NMl 

C 

C  SEARCH  FOR  THE  LARGEST  ENTRY  IN  COLUMN  I,  ROWS  1  THROUGH  H 
C  IPIVOT  IS  THE  ROW  INDBX  OF  THE  LARGBST  ENTRY 
C 

PI VOT  =  0 -  0 
DO  3  J  =  I ,  N 
TEMP  =  ABS { AUG ( J ,  I  ‘  : 

IF< PIVOT. GE. TEMP  GO  TO  3 
P I VOT  =  TEMP 
IPIVOT* J 

3  CONTINUE 

IF  PIVOT. BQ. 0 . 0  GO  TO  13 
I F  <  IPIVOT. EQ.  I  •  GO  TO  5 
C 

C  INTERCHANGE  ROW  I  AND  ROW  IPIVOT 
C 

DO  4  R= I , NP 
TEMP*  AUG (  I,  K  : 

AUG ( I ,  K ) =  AUG  <  IPIVOT. K) 

AUG '  IPIVOT. K) =T8MP 

4  CONTINUE 
C 

C  ZERO  ENTRIES  (  I  *  1 )  ,  (  l  *  2  :  ,  .  .  .  .  ' N ,  1 ,  IN  THE  AUGMENTED  MATRIX 
C 

5  I  P  1  I  *  1 

DO  7  R  -  I  P  1  ,  N 
Q- -AUG  K , I  AUG' 1,1 
AUG  K.  I  0 . 0 
DO  6  J= I  PI , NP 

AUG  K.  J  > =Q*AUG( 1,J j+AUGtl.J) 

6  CONTINUE 

7  CONTINUE 

8  CONTINUE 

if;aug<n,N).bo.o.O)  go  to  13 

c 

C  BACHSOLVE  TO  OBTAIN  A  SOLUTION  TO  AX=B 

C 

X(N'  =  AUG  (  N  ,  NP  )  /AUG  iN.N) 

DO  10  K  =  1  ,  NM 1 
Qr  0 . 0 

DO  9  J- 1 , 8 

Q=Q+AUG(N-K, NP-J)*X(NP-J ! 

9  CONTINUE 

X  N-K  =  •  AUG  (  N-  K  ,  NP  •  -0  AUG  •'  N  -  K  ,  N-  R  1 

10  CONTINUE 
C 

C  CALCULATE.  THE  NORM  OF  THE  RESIDUAL  VECTOR.  B-AX. 

C  SET  IEHROR-1  AND  RETURN 
C 

RSO-0 . 0 

DO  13  I  =  1  .  N 
Q  -  0 . 9 

DO  11  J  -  1  ,  N 


120 


11 


Q*0*A(X.J)*X( J> 

CONTINUE 
BBS  I r  B { I ) -Q 
BMAG  =  ABS  ;  RE S  I 
RSQ=RSQ*RMAG**2 

12  CONT 1  NUB 
RN08M=SQRT( RSQ ) 

IBRROR= 1 

C 

C  ABNORMAL  RETURN  -  REDUCTION  TO  ECHELON  FORM  PRODUCES  A  ZERO 

C  ENTRY  ON  THE  DIAGONAL.  THE  MATRIX  A  MAY  BE  SINGULAR. 

C 

13  I  ERRORS  2 
C 

RETURN 

END 


C 

CM 


mttitut 


CM 


C*  SUBROUTINE  BCM 

CM 

C*  THIS  SUBROUTINE  TARES  THE  KNOWN  BOUNDARY  CONDITIONS  STORED 

CM  BY  THE  PROGRAM  INPUT  IS  THE  FILE  BCM. DAT  AS  WEIL  AS  THE 

CM  " UNKNOWN "  BOUNDARY  CONDITIONS  RETURNED  8Y  THE  SUBROUTINE  GAUSS 

CM  AND  CPEaTBS  THE  8X1  MATRIX  OF  BOUNDARY  CONDITIONS  BCMAT 

CM 


CM 


c 


SUBROUTINE  BCM ( POS , BC , X . BCM AT 
C 

IMPLICIT  REAl*8  (A-H.O-Z'C 
RE  A  L  *  8  8C{4) , X(4> . BCMAT  8) 
INTBGER  POS  <  8 ) 

C 

K  :  0 

L  =  0 

DO  120  1=1,8 

IF  (POS'I)  .GT.  0;  THEN 
K  =  R-*  1 

BCM AT (  n  =X(  R 

ELSE 

l r  L  *  1 

BCMAT  I  r  BC  l. 

ENDI  F 
120  CONTINUE 
RETURN 
END 
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SUBR  UTINE  SOLVER . XMT , BCMAT 


IMPLICIT  RE AL*8  (A-B.O-Z' 


RE  AL*8  GRNMT1  4 , 6  •  ,  GRNMT2  4,6',  BCMAT  fi  .  DT  50,4  ,  SM  5(>  .  4  ,  V 
♦  XMT(50) ,STTAV(50,2  i 

INTBGER  M , N 
$  INCLUDE :  'COMMON. FOR1 
C 
C* 


1 
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c 

Cl 


Ct 

Ct  THIS  LOOP  CALCULATES  THE  LOAD  TERM  A SD  THE  PLASTIC  MO HE  XT 

Ct  TERM  FOR  EACH  EQUATION.  I  REFERS  TO  THE  RADIAL  STATION  AND  J 
Ct  REFERS  TO  THE  EQUATION,  WHERE  THE  FIRST  EQUATION  IS  FOR 
C*  DEFLECTION.  THE  SECOND  FOR  DW  DR.  THE  THIRD  FOR  D"2*  DK~2  AND 
Ct  THE  FOURTH  FOR  D'JW/DR'3-  THE  DERIVATIVES  ARE  THEN  USFD  TO 
Ct  CALCULATE  SLOPES,  MOMENTS  AND  SHEARS 
C» 


DO  155 

C 

IF  ( 1  . EQ.  1  THEN 
RO - B ♦ AD J ( 1 
E  N  o  r  F 

IF  'I  . EQ.  i N ♦ 1 j  :  THEN 
RO  =  A  - AD J (  1  : 

ENDIF 

DO  150  M= 1 . NINT; ADJ , 2 

R  -  B  -DB  L  2 . 0  *■  R  E  A  L  M  tCEL 

CALL  6RNFNC  R . RO , NU . D , GRNMT 1 

STTAV . M, 1 . =  STT  M.9  *STT  M* 1 , 9,  ■  /  2  0 

STTAV  M.2  s; STT' M, 10  -STT  M-1,10  '  2  0 

DO  140  J  -  1  ,  4 

DT  I.J  - DT  I . J  *  XMT  M  *  G  R  S  M  T  1  I  . J  *R  D E L *S TT A V ' M .  1 
-GRNMT 1  3 , J  D-Nt  R*GRSM7 1  2.J-  *R 
-DELtSTTAV  M.2  *G  HNMT 1  2 . J 

1 4  C  CONTINUE 

150  CONTINUE 

IF  I  .EQ.  1  THEN 
RO=B-  DEI. 

ELSE 

RO - RO - D  E  L 
ENDIF 
155  CONT1NLE 
C 
C 

C* 


*  < 


C  * 

Ct  AT  EACH  INCREMENT  OF  WIDTH  DEL  ALONG  THE  PLATE  RADIUS,  CALCULATE 
Ct  THF  DEFLECTION.  SLOPE.  MOMENT  AND  SHEAR.  WHICH  ARE  STORED  IN 
Ct  STT  E  Ar  H  COLUMN  REPRESFNTS  RO .  V.  TH.  M.  U ,  V'1,  V‘‘*  CjfrKES- 

C*  PONDING  T:  EACH  STATION  I  ACROSS  THE.  RADII  S 

C» 

CtitimittmuKtttiiiiiifiMmnmiMitmiifiMtmmnimttii 


C 


DO  200  1  -  1 , N •  1 

IF  I  .  EQ -  ■  THE \ 


123 


B0-8- ADJ (  1  ) 

S NO  IF 

IF  (I  -  80 .  (N-n  )  THEN 
RO :  A  - ADJ I  1  • 

8NDIF 

C 

8  -  A 

CALL  GRNFNC; R, RO , SU, D.GRNMT1  ■ 

R  =  B 

CALI.  GRNFNCl  R,  RO  ,  SU,  D.GRNMT2  . 

00  180  1 . 4 

DO  170  81,4 

SMi  I , J : r AAGRNMTl  .  K , J ; *  BCMA  T  K  -  1 . 0  ** K - 
*  B*GRNMT2  K.J  *8CMAT  K-4  • ' -  1  . 0 J *  * R * SM  I.J 

170  CONTINUE 

SM'  I , J ) =  SM '  I  , J  i  *DT '  1  . J , 

180  CONTINUE 

C 

C  • 


C* 

C* 

C* 

C* 


c* 
c« 
r  9 


STT(Ifl;  -  radius  of  interest 

STT ( I , 2 .  -  DEFLECTION 

STT ( I , 3 )  -  DW  /  DR 

STT; 1.4)  -  RADIAL  MOMENT 

STT: I , 5 •  -  SHEAR 

STTi 1,6)  -  TANGENTIAL  MOMENT 

$TT(I,7;  -  DW'C  DR '2 

STT  1,8-  DW'3.  DR'  3 

STT- 1,9  *  RADIAL  PLASTIC  MOMENT 


STT  1,10-  TASCENTIAL  PLASTIC  MOMENT 


STTi  I  .  1  i  -  HO 

STT.  I  ,2  ’SM:  I  .  1  /RO 

STT-  I.3l=l  0  HO«  SM.1,2  -  STT'  1,2  i. 

STT  . 1.7 )-1.0  RO*SM  1,3  -20ISTT. 1.3  * 

STT  I . 8 :  -  1  .  0  - RO* - SM •  1,4  -3.0«STTiI,7 
STT  1.4:  -  -  D  *  :  STT.  I  ,  7  <  ♦  VU RO*5 TT  .  2.3  ' 

STT  1,5  D*  ■  STT:  I  .  8  *  l  .  0 -RO*STT  '.  I  .  7-  1.0/RC**2.0* 

STT  1.3 

STT  1  ,  6  ■  =  -  D  *  1.0/RO*STT  I , 3 : *  NU*STT :  1,7)  1 


C*  THIS  IS  A  STEP  TO  PREVENT  THE  INTERIOR  POINTS  FROM 

c*  having  the  epsilon  thrown  in. 

c* 

c 

IF  I  - EQ  1  THEN 
RO : B  -DEL 
ELSE 

RO - RO -DEL 
E  N  D  I  F 
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:  0  0  r  L  H  T  :  N  V i 

RETURN 


UMitiiimtutmttti 


THIS  SCBROCT  ZNE  T  A  K  f  S  THE  ELATE  PENDING  SOLUTION  FOR  THi-  ISFf 
LOAD  INCREMENT  AND  CALCULATES  THE  LOAD  At  WHICH  YIELDING  OF  THE 
OTTER  F'BEfi  *  I  L  L  OCCUR.  THE  LOAD  MATRIX  IS  THEN  INCREASED  TC 
THIS  VALLE-  THIS  IS  ESSENTIALLY  A  TIMF  SAVER.  TO  AVOID  STEFFIS 
IT  THE  LOAD  IN  THF  ELASTIC  RANGE . 


urn 


its 


itiiiiniiiifnttioitiimiiiimMMttttimiit 


i  .  P  HOI  TINE  Y 1  I>  KK . 1MAX . LDMaTI 

implicit  k  e  a  i.  *  y  a-h.o 

HML‘1  ID^aTI  so  . STRMaX . S’TMAX 
:  NTEGLn  KR  .  I  MAX  .  I  HM.A  X  .  I  TMAX 

I  ■  r-F  'I*  •  M M N  F'jR  ' 


FIND  THI-  LOCATION  ACROSS  THE  PLATE  WHERE  THE  BENDING  MOMENTS 
ARE  THE  HIGHEST 


STRMAX  0.0 
S  .-tmax  0  0 

I  ■’  k  2L  I  1  .  SI  NT  A D J  2  •  I 

IF  ABS  STT  1.4  G7.  STRMAX  THEN 

STRMAX  APS  STT  1.4 

: rma.x  •  ; 
f  n  ,  :  i 

;  i  a  .**  ' '  : .  ►  :  .  :  :  *  a  »  :  h  1  n 

•  ■  '  ah  - : :  : . 


■•'►“o  i :  ? : : *  o  then 

A  F  . 
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lt  =  2 

imax=:tmax 
END  I  F 

122  CONTINUE 

C 

C* 

C*  CALCULATE  THE  EQUIVALENT  STRESS  AT  THE  LOCATION  WHERE  STRESSES 
C*  ARE  HIGHBST.  COMPARE  IT  TO  THE  YIELD  STRESS.  THEN  RATIO  IP  TKL 
Ct  LOAD  MATRIX  ACCORDINGLY 

C* 

c 

ERMAX  =  STT ;  IMAX.T  tGEOMtS  .'2.0 

BTMAX=STT ' IMAX , 3 • *G£OM  5 .  '2.0  ST7  I MAX, 2 

SRMAX  =  GEOM  6  i  -GEOM'  7  .  .  1  .  O-GEOM  3  **2.0  »  ERMAX • GE  OM  3  *E?WAX 
STMAX=GEOM ' 6 WGEOM  7  '1  O-GEOM  3  **20  *  ETMAX-GEOM  3  *PRMAX 
SEMAX  =  SORT  SRMAX**2. 0- S RMA X *S TMA X  - S TM A X *  * 2  -  0 
R -GEOM : 6  •  /SBMAX 
DO  126  J= 1 . NINT' ADJi 2  • 

LDMATI  J  ■ =  R*  LDMA T I :  J 
226  CONTINUE 

WRITE; * . 130 )  LDMATI i 1 • 

130  F0RMAT(2X. ’ LDMATI  =  ‘.F20.I0 

C 

RETURN 

END 


C 

C* 


*  SUBROUTINE  PLAST 


v.  v 

C* 

C 


THIS  SUBROUTINE  CALCULATES  PLASTIC  MOMENTS  USING  A  TRAPAZOIDAL 
RULE  INTEGRATION  SCHEME.  THE  SUBROUTINE  WILL  ONLY  HANDLE 
LINEAR  ELASTIC-PBRFECTLY  PLASTIC  MATERIALS. 


*********************************************** 


SUBROUTINE  PLAST'  1 


C 

IMPLICIT  HEAL*P  AH. 0-7  C 
RE  A 1*8  MRP, M TP. Nl 
INTEGER  IZ 

S INCLUDE  ’ COMMON •  FOR  1 
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f 


c 

ct 
ct 

c*  variables 

c* 

C*  DER 

C*  DBT 

C*  DERE  “ 

Ct  DETE  - 

C*  DEL 

Ct  I 

C*  IZ 

Ct  DSl 

Ct  DS2 

C*  STRT  - 

Ct 

Ct  SET 

C* 

Ct  SE 

Ct  DEPT  - 

Ct 
C* 

ct  ft 

ct 

ct 

ct  z 

Ct 

ct 

c 

STT-  I , 9  .  =  0. 0 
STT ; 1,10=00 
SSUMR -  0 . 0 
SSt'MT  =  0. 0 

DEL=GEOM  < 5 ; / 2 . 0 / ( ADJ ' 6 ) -  1 . 0 • 

Z  =  0 

E=GE0M<6) /GE0M(7) 

NU  =  GEOM  y  3 ? 

C 

DO  500  IZ  =  NINT: ADJ ( 6  ‘  ’  , 1 , -  1 
C 

c* 

c* 

Ct  THIS  SECTION  CALCULATES  THE  STRAIN  AND  STRESS  INCREMENTS 
Ct  CORRESPONDING  TO  THE  GIVEN  LOAD  INCREMENT.  IT  THEN  ADDS  THE 
Ct  STRESS  INCREMENT  TO  THE  PREVIOUS  TOTAL  STRESS  TC  CALCULATE  THE 
«"•  ECL  I  VALENT  STRESS.  THE  EQUIVALENT  STRFSS  IS  COMPARED  TO  THE 
Ct  UNIAXIAL  YIELD  STRESS  TO  SEE  TE  YIELDING  HAS  OCCURRED.  IF  SO. 
Ct  THE  PLASTIC  STRAIN  INCREMENT  IS  SET  EQUAL  TO  THE  TOTAL  STRAIN 
'*  INCREMENT 
C  t 


USED  IN  THIS  SECTION  ARE  DEFINED  AS  FOLLOWS 

TOTAL  RADIAL  STRAIN  INCREMENT 
TOTAL  TANGENTIAL  STRAIN  INCREMENT 
ELASTIC  RADIAL  STRAIN  INCREMENT 
ELASTIC  TANGENTIAL  STRAIN  INCREMENT 
DBPTB  INCREMENT 

RADIAL  POSITION  STATION  NUMBER 
HALF- THICKNESS  STATION  NUMBER 
RADIAL  STRESS  INCREMENT 
TANGENTIAL  STRESS  INCREMENT 

MATR I X  HOLDING  THE  TEMPORARY  STRESSES  AT  THE 

GIVEN  LOADING  INCREMENT  ■;  l  =  RADlAi.  2  =  TANGENTIAL  . 
TEMPORARY  EQUIVALENT  S  TRESS  -  USED  FOR  YIELD 
CRITERIA 

EQUIVALENT  STRESS  FROM  PREVIOUS  INCREMENT 
ARRAY  CONTAINING  THE  TEMPORARY  STRAIN  INCREMENTS 
AT  THE  GIVBN  LOADING  INCREMENT  U=RADIAL;2= 
TANGENTIAL 

CONSTANT  REFERRING  TO  THE  FRACTION  BEYOND  YIELDING 
THAT  OCCURRED  DURING  THE  COURSE  OF  THE  LAST  LOAD 
INCREMENT 
DEPTH 
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c 

DER- -STT : I . 7 ' •  I Z-  1  *DEL 
DET--STT; 1,3  1Z-1'*DBL'ST{  1,1  • 

DERE  =  DER -  DEPT ■  I.IZ.l, 

DE  TB  =  DBT-DEPT .  I,IZ.2‘ 

DS1  =8  /  (  1 . 0  -  N  U  *  *  2 . 0  ,  *  •  DERE-NL’*DETE 
0 S 2 a E /  (  1  ■  0-Wt?**2  0  *  ,  DETB*  Nt’*D8RE 

strt:i,iz,i:=str. i , iz, l ;*dsi 

STRT  :  I  ,  IZ  .  2  -  =STR  I,2Z,2>DS2 

SETi  1  ,  IZ  1  -SORT  STRT  .  I.IZ.l  »*2 . 0-STRT.  I  .  IZ,  1  ‘  *STRT  I  ,  IZ,  2 
*  ♦STRT:  I  ,  IZ,2‘,  **2  . 

c 

c* 

C«  SEE  :f  yield  criteria  is  exceeded,  if  so,  rezero  the  plastic 

C*  STRAIN  INCREMENTS  AND  THE  STRESS  INCREMENTS  ALSO,  RETURN  TO  THE 
C*  MAIN  PROGRAM  IF  THE  OUTER  FIBER  HAS  NOT  BEGIN  TO  YIELD 

C* 

C 

IF  v'  SET  '  I  ,  IZ  ;  .IT.  GEOM  B)'  THEN 
DBPT(  I  ,  IZ,  1  )=0 . 0 
DBPTf I ,  IZ,2  *  =0. 0 
STT { I , 9 ) =0 . 0 
STT I  I , 1 0 ' =  0 . 0 

IF  (IZ  .80.  NINT<ADJ!6M)  THEN 
RBTURN 
ENDIF 
ELSE 

SET:  I  ,  IZ  -GEOM  6 
DEPT  I  ,  IZ,  1  • =DER 
DEPT;  I  ,  IZ . 2  •  -DET 

IF  i SE :  I . IZ  . LT.  GEOM  6  »  1  THEN 

R=(SET  I . IZ ■ -GEOM  6 . : / : SET! I . IZ  -SE'I.IZ  : 

DEPT  I.  IZ,  1  >  *  R  *  D  £  P  T  I.IZ.l 
DEPT: I, IZ . 2  ; R*DEPT  I, IZ,2 
STRT  I.IZ.l  :=STRi  I.IZ.l  '  *  DS  2  *  I-R; 

STRT ( I.I2.2-STR. 1,IZ,2'*DS2»: 1  -  R  1 
ENDIF 
ENDIF 
C 

Z=  Z ♦ 0 E  l 
C 

500  CONTINUE 
C 

C* 


PERFORM  A  TRAPAZOIDAl.  RILE  INTEGRATION  OF  THE  STRAIN  OVER  THE 
PLATE  HALF  THICKNESS .  AND  MULTIPLY  BY  2  TO  GIVE  THE  TOTAL  RAIIaI 
AND  TANGENTIAL  PLASTIC  MOMENT  INCREMENTS  STT  I .9  AND  STT  J . I f 


C* 

C 

Z - -DEL  2-0 

DO  54  C  1 2  ^  1 . NIST  ADJ  6  -  1 
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Z  =  Z  *DE  L 

DMPR=  <  (DEPT<  1 ,  IZ.  1  •  *DSPT(  I ,  IZ-*  2  ,  1  j  2-GEOM:  3  ;  •  .  DEPT  '  I .  IZ ,  2  , 

*  *DBPT<  I,  IZ-1.2:  >J*Z/2.0 

DHPT=C , DBPT{ I . IZ,2)*DEPT< 1 , IZ*1 ,2; )*GEOM  3  : DEPT, I , IZ , 1 

*  *  D  EPT ( I . IZ» 1 , 1 ) )  , «  Z / 2 . 0 
SSUMR  =  SSUMR  *•  DMP  R 

SSUMT  =  SSUMT ♦ DHPT 
540  CONTINUE 
C 

STT.  I , 9; =2 . 0*GEOM' 6) /GEOM' 7)  /  { 1-GBOMr  3  **2.0*DEL*SSIMR 
STT( I , 10 ; =2. 0*GEOM( 6  I  /GEOM  (  7  )/(l-GBOMt3)**2.0  / *DEL*SSUM7 

WR I TE i * , GOO  '  I 

600  FORMAT  2X LAT  POS  *  ’.12 

WRITE*:*, 700  STT.1,9  ,STT  1,10 

700  FORMAT' 2X. ‘MRP  =  * ,F20.12,‘  MTP  =  ’ .F20.12 
C 

RETURN 

END 


C 

C* 


c* 

C*  SUBROUTINE  UPDATE 

C* 


C*  THIS  SUBRCU71SE  TAKES  TBS  INCfiEMEN  TAl  STRESSES,  STRAINS, 

DEFLECTIONS.  MOMENTS.  ETC  AND  ADDS  THEM  TO  THE  PREVIOUS  TOTALS 
ONCE  CONVERGENCE  OF  PLASTIC  MOMENTS  HAS  OCCURRED. 


C* 

C 

SUBROUTINE  UPDATE l STT 1 . IMAX , OUT . LDMAT 1 . LL. LDMATO 
C 


IMPLICIT  RE AL*8  (A-H.O-Z 

RE AL*8  STT 1(50, 2) , LDMATO ; 50 ), LDMAT 1 < 50 , . Nl , DEL.E 
INTEGER  NN , OUT , LL 
(INCLUDE . 'COMMON. FOR' 

C 


NN-NINT  ADJ  2  *1 

E - GEOM  6 •  GEOM  7 
DFL-GEOM  5  20  ADJ  6  i  0 

- GEOM  3 


DO  i 000  I  - i  . NN 
DO  905  J: LL.  10 

ST  I.J  -  S  T  I  ,  J  -STT  1  .  J 
905  CONTINUE 
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LDHATO { I i-lDHATi . I  •♦IDMaTO,  I  . 

C 

DO  980  12=  1 , N I  NT i A  0  J  6 

Ci 

C*  If  THE  YIELD  CRITERIA  BAS  BEEN  EXCEEDED.  THEN  ADD  IN  THf 
C*  PLASTIC  STRAIN  INCREMENTS  AND  UPDATE  THE  TOTAL  STRESS  TO 
C«  EQUAL  THE  TEMPORARY  STRESS  CALCULATED  IN  THE  SUBROUTINE 
C*  PLAST 

C* 

C 

IF  ,  S  B  '  l  ,  12)  -GT.  G B OM ' 6 :  >  THEN 

8  P  I ,  IZ.  1  ‘  =EP-  I ,  IZ ,  1  ' *DBPT  I.IZ.l 
BP'  I.  IZ , 2 ,  ‘ EP i  I,  IZ.2  -DEPT  I.  12.2 
STR  I ,  IZ .  I  ■  =  STRT -  I.IZ.l 
STR,  I , IZ.2 • =STRT  I  .17 .2 
Z= < CBOM; 5 > /2- 0/ I ADJ< 6  -  1 . 0  t  17-1 
ELSE 
C 

C* 

C*  IF  THE  YIELD  CRITERIA  HAS  NOT  BEEN  EXCEEDED.  THEN  UPDATE  THE 
C«  STRESSES  AND  EQUIVALENT  STRBSSBS  USING  THE  ELASTIC  EQUATIONS 

C* 

C 

STR  I  ,  IZ.  1  )  =  -E/ (  1 . 0-NU**Z . 0 , *  ST  I . ?  t  l 7- i  •  *0E1 
8  P  (  I,  IZ,  1  S  ♦  N  U  *  '  S  T  >'  1,3  *  IZ-1  )*DEL.  ST,  I  ,  1 
STR  I , IZ, 2) --t/ £  I . o- NU**2. 0 , * : S T  l . 3  j  » .  I Z- I  •  DEL  ST 
•  8  P  <  I, IZ,2»*NU*:ST. 1,7 )•;  12  —  1  ;  *  D  E  L  ’ 

SB.'  I  .  IZ  -*SQRT  STR  I,  IZ,  J  <  **2 . 0- STR  <  I.  IZ.  I  :* 

STR  I  ,  IZ.2;  ‘STR:  ! .  IZ.2  *#2 

END  I  F 

IE  v  ABS  EP;  I  ,  IZ.  1  /  •  .GT.  GEOM  9  •  OR.  ABS  EP  I.IZ.L* 
.GT.  GEOM  9  THEN 
OUT=  1 
E  N  D  I  F 

980  CONTINUE 

1000  CONTINUE 
C 

LL=  2 

RETURN 

BND 
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c 

c* 

c* 

C*  SUBROUTINE  RBSUL 

C* 

C*  STORE  THE  H8SULTS  AT  THB  8ND  OF  EVERY  10  LOAD  INCREMENTS  IN 

C*  THB  FI  LB  OUT. DAT.  ALSO  SCREEN  DUMP  THBSB  RESULTS. 

Ct 

Cttttttttttt*tt**»**t9****t**ttttttt#*t*sttttttt**ttt**t»t*t*ttt*t**» 

Ct 

c 

SUBROUTINE  RBSUL ( MM , STT 1 , COUNT , LDMATO ) 

C 

IMPLICIT  RE ALt8  (A-H.O-Z) 

RBALt8  STT1 (50, 2) , LDMATO(50) 

INTBGER  NN 

t  INCLUDE :  'COMMON. FOR' 

C 

NN=NINT ( ADJ ( 2 ) )  +  l 
DEL=(GEOM( 1)-G80M(2) )/ADJ(2) 

C 

IF  (COUNT  .BQ.  1.0)  THBN 
WRITBf 15, 100) 

WRIT8(t, 100) 

100  FORMAT (  /  //  ,  ’**i**»»»**t»im**l»***mM*****t*****M*in*i» 

+  ’***$**««** **t**t**t*t*t*«**«**t*‘,//) 

WRITE( 15, 120)  MM, LDMAT0( 1 )/DEL 
WR I TB ( t , 120)  MM, LDMATO ( 1 )/D£L 
120  F0RMAT(2X,  *  LOAD  INCH  =  ' , 1 3 ,  *  LOAD  =  ’,1,0.4,/) 

C 

WR I TE ( 15, 130) 

WRITE ( * ,  130  ) 

130  FORMAT(2X, ’ RADPOS ', 2X, ’ DEFLECT  ’,2X, ’SLOPE  ’,2X,’RAD  MOM  ’ 

♦  , 2X ,  ’  MRP  * , 2  X ,  ’  TAN  MOM  ’,2X,'  MTP  SHEAR  ’,/) 

DO  160  1=1, NN 

STMl =ST( I , 4 ) *ST( 1,9) 

STM2=ST( I , 6 ) +ST( I , 10) 

WRITB{ 15, 140)  ST( 1 , 1 ) , -ST( 1 ,2) ,ST(  I , 3 ) , STMl , ST<  I , 9)  , 

♦  STM2,ST( I ,  10) ,ST( 1,5) 

WRITB( t, 140)  ST ( 1 , 1 ) , -ST ( I.2),ST(I,3),STM1,ST(I,9). 

♦  STM2 , ST ( I ,  10) ,ST( 1,5) 

140  FORMAT (2X , F5.2,2XfF7.5,2X,F7.4,2X,FlO.O,2X,F5.O,2X,FlO.0, 

♦  2X,F5.0,2X,F10.0) 

160  CONTINUB 

C 

WRITE ( 25, 163) 

WR I TB ( t , 163) 

163  FORMAT ( / , 2X , 

♦  ’ _ i;... 

DO  1000  1=1, NN 

IF  ( ST ( 1,9)  .NF.  0.0  .OR.  ST ( I , 1 0 )  . NE .  0.0)  THEN 
WR I TE (  15, 170)  ST( I. 1) 

WR  I  TB  (*,170)  STU.l) 

170  FORMAT(/ ,2X,  ’ RADIAL  POSITION  =  ’  , F5 . 2  ) 
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} 


165 


180 


200 

1000 


BRP 
•  ,ZX,  ' 


.  2X  , 
S  E 


SilAD 

*  ,/  > 


WRITE { 15,165) 

WR I TB ( * , 165)  , 

PORMAT ( / , 2X , ’  DBPTB  ’  ,  2X  , 

2x, •  8TP  ' , 2X , '  STAN 
DO  200  IZ=1.NINT(ADJ(6) ) 

IF  (SB(I.IZ)  -GT-  GBOM ( 6 ) )  THEN 

Z; ( GBOM (5)/2.0/(ADJ(6)-1.0))*(l2-ly 

WRITS'  15, 180)  2fBP(l,IZ,l).STR(l.  12.  D.BPd,  IZ.2). 

STR{1.1Z,2),SEU,IZ)  7 

WRITB(*.  180)  2.FPU.  IZf  1)  ,STR(I.  I*.  l>  ,BP(I.1Z.2J  • 

STR(1. 1Z.2)  ,SBU  ,  12)  o  rfi  ! 

FORHAT(2X,F6.4.3X,n0.8,3X,F8.  1 , 3X  ,  FI  0 . 8 . 3X  .  T8  .  1. 

3X , F8 ■ 1 ) 

bndif 

CONTINUE 
END  I  F 
CONTINUE 
BNDIF 
RETURN 
END 


SSS  *0>.O.PTC50.»0.2t.S«<50.50.= 
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